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SOLUTION OF ELASTIC AND ELASTO-PLASTIC PROBLEMS 
BY THE METHOD OF LINES 


Abstract 

by 

ALEXANDER MENDELSON and JAVED ALAM 

An improved formulation of the method of lines (MOL) Is pre- 
sented. The line method lies ml<: : v between completely analytical 
methods and completely discrete methods such as finite differences. 
The five-point finite difference nnulas are Introduced to reduce a 
set of coupled partial differential equations Into a set of simul- 
taneous ordinary differential equations. The resulting ordinary 
differential equations are solved by a recurrence relation method, 
which is found to be very suitable in solving two-point boundary 
value problsims . 

The formulation is further extended to include the small scale 
plasticity effect. The Von-Mises criterion is used for yielding 
and isotropic hardening rule is followed to decide for subsequent 
yielding. Prantl-Reuss equations are employed as constitutive re- 
lations, The resulting nonlinear equations are solved by using the 
method known as successive elastic solutions. 

For two specific geometries namely the edge notch specimen and 
the compact tension specimen, the complete field solutions for 
stresses and strains are obtained. 
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A numeirical experimentation is carried out to establish the 
convergence characteristics of the improved MOL. The results for 
stress intensity factors (SIF) are also compared with the existing 
results obtained by finite element analysis. 

The stress intensity factors are calculated for the compact 
tension specimen containing curved crack fronts. A thickness aver- 
age stress intem-ity factor is evaluated and it is used to compare 
results for fracture toughness. 

To explain the tunnelling behaviour in the compact tension 
specimen a complete elasto-plastic analysis is carried out, using 
the improved formulation. For a straight crack the load versus 
crack mouth opening displacement curve obtained by analysis is com- 
pared with the experimental plot. The experiments were conducted 
at NASA Lewis Research Center. The J-integral values at different 
thickness levels of the compact tension specimen are computed to 
predict the crack Initiation through the thickness. 

The results obtained show that the introduction of the five- 
pclnt finite difference formulas has considerably enhanced the ac- 
curacy of the method of lines. A relatively coarse grid is suffici- 
ent to yield an accurate result. 
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CHAPTER 1 
INTRODUCTION 

Knowledge of the stress and strain distributions in the neigh- 
borhood of a slnguj^arity such as the tip of a crack In a beam or 
compact tension specimen Is of fundamental Importance in evaluating 
the resistance to fracture of structural materials. In the past 
various researchers (lt2) using stress functions and complex variable 
approach, solved plane elastic problem for cracked bodies. Gross 
et. al. (3) used boundary collocation technique to solve different 
planar crack problems. The success of both analytical and numerical 
methods^^have placed present fracture mechanics on a firm foundation. 

However, many fracture mechanics applications do not Involve 

the solution of plane problems. For example, the stress analysis of 

corner cracks or surface cracks In practical structures is a tnree- 

dlmenslonal problem and rational approach to design of such structures 

requires an accurate method of three-dlmeqslonal analyses. As de- 

% 

scribed In Reference A, at present, there Is a wide divergence among 
the various approximate solutions to such three dimensional crack 
problems and in these circumstances structural design using fracture 
mechanics must proceed largely on the basis of experience and labora- 
tory modeling tests. 

A closely related problem occurs during the fracture toughness 
testing of standard compact tension specimens as shown In figure 1(a). 
Very rigid standard have been Imposed by ASTM (5), to ensure that 
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ViiHvI fi'AOtur« tv>u|khn««it vdlut#* will b» obt«ln««i In a glvan taat. On* 
of tha raqulramanta ia tha production of a fatltua crack. Howavar. 
tha fattgua crack will uaually not «row uniforvly acroaa tha apaclman 
thicknaaa» l.a. tha Initially atraight crack front will bacoiaa curvad. 
Tha crack lanitth will thua vary acroaa tha apacinan thlcknaaa, and 
aoma avaraga valua nviat ba uaad to calculata tha fractvira toughnaaa. 
Kataranca provtdaa a atandard waaaura auch that if tha crack 
front curvatura aa taaaavirad afiar tha apaciman ia brokan, ia graatar 
than thia atandard awaanra, tha taat ia Invalid. 

Kacantlv. howavar. auygaativ'na hava baan niada (b) that tha 
abova atandard la too rigid and can ba ralaxad ainca tha affact of 
tha curvad crack front may ba laaa than origlnallv anticlpatad. Thaaa 
auggaatlona, howavar, at a baaad on rathar waak avtdanca t?) and it 
iii not vat claar whathar tha K- I*!*! atandard ahould ba ralaxad or not 
In thia raapact . » at naadad la a claar daflnltlva analyala of 

tha at tact of tha crack front curvatura on tha atraaa intanalty fac- 
tor. ihta of tha objactlvaa of tha praaant invaat igat ion ia to par- 
form thia ana I va la. 

Initial attampta to aolva fvillv thraa-dimanaional alaatlc 
proMaam wara mada bv applying cartaln bovutdary corractlon fact ora 
to tha axiating tw\> dinanaional raaulta. Thaaa approximat Iona wara 
aolav baaad on intuit iva Judgamant and raaultad In wldaly acat - 
tarad aolnt Iona tor tha aama problama. Thraa-dimanaional problama 
aia inhalant Iv tix'ra complax and in aloK'at all tha caaaa ona haa to 
•aaort to aomc sort of nuraarlcal tachx'iquc which laada tv' a largo 
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set of equations to be solved on a digital computer. Many numerical 
schemes are devised to obtain the stress and strain distribution in 
the neighborhood cf a crack in a three-dimensional body. Notably 
among them are Boundary Element Method (8) and Finite Element Method 
(9,10,11). Significant progress has been reported in the develop- 
ment of both of these techniques. Special crack tip elements are 
developed which include the effect of crack tip singularity. However, 
the effect of element size and arrangement around the crack tip, on 
the convergence of the solution needs more extensive study. 

A few attempts were made to determine stress fields near 
crack tips having some curvature by using experimental techniques . 
McGowan (12) used a 3-D stress freezing technique to investigate 
the effect of crack front curvature on the stress intensity factor 
distribution in a single edge notch specimen. He reported that the 
increase in the crack front curvature decreases the stress intensity 
factor (SIF) at the center of the specimen and also the thickness 
average stress intensity factor. In his study the value of SIF at 
the surface was insensitive to the change in crack front curvature 
as well as the Poisson’s ratio. This study predicts the behavior, 
which is In good agreement with other analyses, but the actual results 
have discrepancy of about 35% with an existing finite element analysis. 
The discrepancy in the results has been attributed to the difference 
in the crack length in the two analyses. 

Fourney (13) devised a new technique of analyzing three-dimen- 
sional problems utilizing scattered light speckle Interferometry. He 
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also found that the stress intensity factor increases at the surface 
due to the presence of crack front curvature. The crack tip singu- 
larity was found to be strongest at the center. This was cited as 
the reason for tunnelling behavior in a compact tension specimen. 

Peirera et. al. (14) carried out a finite element stress analysis 
of a compact tension specimen. They distorted the crack tip elements 
to model a circular type of crack front. It was found that the local 
values of SIF decreases at the center and it rises at the surface of 
a compact tension specimen. An average value of SIF based on energy 
considerations was also calculated. This value was found to be de- 
creasing with increasing crack tunnel depth. 

Neale (15) performed another finite element analysis for a non 
standard compact tension specimen (W/B « 8) containing a thumbnail 
crack. He showed that if the average crack length is calculated 
based on ASTM method, then the fraccure toughness is overestimated. 
Both the previous analytical studies did not use the standard compact 
tension specimen as shown In figure 1(a). Peirera et. al. used a 
tensile loading which does not model the actual loading condition. 

In their studies a crack length to width ratio (aA/) of 0.25 was used, 
while ASTM standard prescribes a range of 0.45 to 0.55. In Neale's 
analysis the thickness of the specimen used was k of the standard 
thickness specified in the ASTM standard. Consequently, the results of 
these analyses can not be directly applied to a standard compact 
tension specimen. 

All these analyse^ predict that crack growth will initiate from 
the surface of the specimen due to the presence of maximum SIF there. 
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This finding, is in direct contradiction with the experimental results, 
in which a tunnelling behavior is observed. Neale (16) tried to ex- 
plain tnis anamoly on the basis of plasticity effects. He performed 
an approximate elasto-plastic analysis to improve the numerical 
values of SIF. His elastic-plastic thumbnail model indicates that 
the maximum stress intensity factor occurs in the central portion of 
the compact tension specimen which is in agreement with experimental 
observations. However, due to the approximate nature of the analysis, 
there is a need to perform an accurate elasto-plastic stress analysis 
to resolve this discrepancy. This is another major objective of the 
present investigation. 

The recently developed method known in the literature as the 
Method of Lines (MOL) was selected to perform the elastic and elasto- 
plastic stress analysis of two specimens namely the edge notch speci- 
men shown in figure 2 and the standard compact tension of figure 1(a). 
In this method there is no prior assumption on the stress field 
around the crack tip. There is an extensive literature, primarily 
of Russian origin for the method of lines, which has been summarized 
to 1965 by Liskovetes (17). The method has been used extensively in 
solving problems in the area of fluid mechanics. 

Jones, South and Klunker (18) presented an analysis of conver- 
gence and stability for the case of linear elliptic partial differen- 
tial equations solved by the MOL. They ind'*cated that the results 
obtained are of sufficiently high accuracy. It was shown that if 
the region of interest is divided into sufficiently few strips by 
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the dividing lines then accurate solutions can be obtained by using 
h.f^her order finite difference approximation. In the study by Kurtz 
et. al. (19) It was observed that In problems where standard techni- 
ques failed to converge the line method was able to produce results. 
Hopkins and Walt (20) compared the execution time of MOL with 
Galerkln and collocation techniques, which were viewed as finite 
element discretization. With very few exceptions, they found that 
the line method yields faster results as compared to the other two 
techniques for a set of coupled and \mcoupled parabolic partial dif- 
ferential equations. 

The line method lies midway between completely analytical 
methods and completely discrete methods such as finite differences. 
The basis of this technique Is the substitution of finite differences 
for the derivatives with respect to all the Independent variables 
except one for which the derivatives are retained. This approach 
replaces a given partial differential equation with a system of 
simultaneous ordinary differential equations whose solution can then 
be obtained by standard means. These equations describe the depen- 
dent variable along lines which are parallel to the coordinate In 
whose direction the derivatives were retained. 

An Inherent advantage of the line method over other numerical 
methods Is that good results are obtained from the use of a relative- 
ly coarse grid. This use of a coarse grid Is permissible because 
parts of the solutions may be obtained in terms of continuous func- 
tions. Additional accuracy in normal stre.-^s distributions is derived 
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from the fact that they are expressed as first-order derivatives of 
the displacements and these derivatives can be analytically evalu- 
ated. Inherently inaccurate numerical dif ferentation is required 
only for evaluating the shear stresses, but this presents no impor- 
tant loss of accuracy in this study since they are usually an order of 
magnitude smaller than the normal stresses. For problems with geometric 
singularities, additional accuracy is derived from using a displacement 
formulation since the resulting deformations are not singular. 

The method has been shown to be well suited to the solution of 
certain three-dimensional crack problems, as has been successfully 
demonstrated (21,22). Fu and Malik (23,24) are the first to report 
the application of MOL to elasco-plasticity . The formulation was 
presented in terms of displacements and their normal derivatives. 

These equations were solved by a combination of power series and 
modal matrix method. 

Although the MOL has given very good results in a number of 
specific geometries as described in the previous references, its 
use has been limited heretofore to bodies with rectangular boundaries 
(including the crack boundaries), or to the problems of axial symmetry. 
Thus, the very important surface flaw problem or curved crack front 
problem have not. been treated. This was primarily due to the limi- 
tation in the tiumber of lines that could be used, since the rate of 
convergence of the infinite series solutions of the differential 
equations decreases as the number of lines increased and as the 
dimensions of the specimen increased. It has also been shown in 
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reference (18) that the use of the MOL for solving elliptic partial 
differential equations can lead to problems of Instability if the 
line spacing Is made too small. 

A third major objective of this Investigation was therefore to 
reformulate the method to avoid or minimize the above difficulties. 
This was done by employing 3- point finite difference formulas rather 
than the usual 3- point formulas and by devising an algorithm for 
solving the three sets of ordinary differential equations for the 
two point boundary value problem using recurrence relations. These 
Innovations have appreciably increased the applicability and accuracy 
of the MOL. The elastic formulation Is then modified by Including 
the plastic strain terms to solve an elasto-plastlc problem. 

A complete formulation In a rectangular cartesian coordinate 
system Is presented for the newly Improved MOL. The governing field 
equations are augmented to Include the plastic strain terms. To 
study the convergence characteristics of the method, complete elastic 
solutions are obtained for different numbers of lines for an edge 
notch specimen. The stress Intensity factors are computed for an 
edge notch specimen and a standard compact tension specimen. They 
are compared with the existing finite element analyses results to 
assess the validity of the MOL. The complete field solutions are 
evaluated for both the specimens for the different cases of crack 
front curvature. For compact tension specimen of a given material 
(A 2,-5083) complete elastic-plastic stress analyses are carried out. 
Similar analyses are also done for specimens with different curved 
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crack fronts. Various fracture mechanics parameter such as the 
Jo Integral are evaluated to establish the effect of crack front 


curvature . 
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CHAPTER 2 

SOLUTION OF CRACK PROBLEM BY THE 
METHOD OF LINES 

In this chapter a method of solution for three-dimensional 
elasto-statlc problems Is described. All the underlying assumptions 
and the governing field equations are described. An Improved method 
of lines, in which five point finite difference equations are used, 
is employed to obtain a set of second order ordinary differential 
equations from the governing coupled partial differential equations. 
The second order differential equations are rearranged to formulate 
a set of first order differential equations and these equations are 
then solved by a new algorithm involving recurrence relations, 
which was found to be a simpler and an efficient technique to solve 
two point boundary value problems. The last section of the chapter 
deals with the evaluation of boundary vectors used in the recurrence 
relations for an edge notch specimen. The details are also given 
for the loading idealization for a compact tension specimen and its 
incorporation into the equations of the method of lines. 

For all the elasto-statlc problems discussed in this thesis the 
following assumptions shall apply. 


a) The deformations are infinitesimal. 

b) All deformations are elastic. 

c) Materials are isotropic and homogeneous. 
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d) Body forces are neglected. 

2.1 Governing Equations 

Within the framework of linearize-^ elasticity theory, the field 
equations neglecting body forces are listed below. 

The equilibrium equations are, 


0 • 0 i,j ■ 1,2,3 

f J 


( 2 . 1 ) 


where the standard tensor subscript notation is used. Hooke's law 


is, 


0 


ij 



e 

kk 


+ 2 G e 




( 2 . 2 ) 


where X and G are Lame's elastic constants. The strain-displacement 
relations are. 


c 


ij 



"j.i> 


(2.3) 


The solution must satisfy these equations at all interior points of 
the body and, in addition the prescribed conditions must be met on 
the boundaries. The above three sets of equations are combined to 
form three partial differential equations in terms of displacements. 
The resulting equations which are known as the Navier equations of 
equilibrium are. 


G Uf + (X + G) u^^j^ - 0 i,j - 1,2,3, (2. A) 

Vor problems formulated in rectangular cartesian coordinates, the 
equations can be written as. 
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- 9 /9v 3w. 

dx By Bt 


7 2 2 

9^v ,3 V . 3 V. 

— 2 + C^ (-^ + -^) 

ay 9z dx 


3 .3u 3ws 

*=2 55 's; + 5J> 


3^ . ,3^ . 3^, 

57 ^ V ^ 


<=2 51 ' 5 ; * 55 ^ 


r . 1 ^ 

1 2(1 - V) 


and ^2 - - 2(1 -^) 


V being Poisson's ratio. 

The stress-displacement relations in cartesian coordinates are. 


®x " (1 + V) (1 - iv) 


r . . /9v , 3w. , 

((1 . „) _+ v(jj + 5i>’ 


V (1 + m) (1 - TvT 


r V ov . /9w . 3uv . 

"vT 55'^'’ '55 ■^55” 


* (1 + v) (1 - 2v] 


55^'' '55*55>’ 


E r 9v . 9u, 

°xy “ mr^ 5^1 


E ,3v . 3w, 
r-TTTT 


2(1 + V) ^3r 5y 


_ _ E r9w . 9ui 

°« JTrrJT ‘55 * 3i' 


(2.U) 
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The method of lines Is used to solve the three coupled partial 
differential equations given by equations (2.5) through (2.7). For 
each problem appropriate boundary conditions are satisfied. In the 
next section the line method is described in detail. 

2.2 Method of Lines 

The recently developed Method of Lines (MOL) has been shown to 
be well suited to the solution of certain three-dimensional crack 
problems. In this method there is no prior assumption on the stress 
field around the crack tip. The line method lies midway between 
completely analytical and grid methods. The basis .'f the method is 
substitution of finite differences for the derivatives with respect 
to all the independant variables except one, with respect to which 
the derivatives are retained. This approach replaces a given partial 
differential equation by a system of simultaneous ordinary differen- 
tial equation whose solution can then be obtained by standard 
means. The equations describe the dependent variable along lines 
which are parallel to the coordinate in whose direction the detiv.-- 
tives were retriined. 

Since in three-dimensional elasticity problems solutions of 
three partial differential equations are desired, three sets of 
parallel lines are constructed. An arbitrary grid consisting of 
these three sets of parallel lines is shewn in figure 3. The lines 

parallel to the x axis are numbered as 1,2,3, — hT, NY •¥ 1 2NY, 

2NY +1 3NY , 3NY + 1 — i. The lines parallel to the y axis are 
numbered 1,2,3, — NZ, i:Z + 1 — 2NZ, 2NZ + 1 — 3NZ, 3NZ + 1 — m. 
Finally lines parallel to the z axis are numbered as 1,2,3, — NX, 



lA 


NX + 1 — 2NX, 2NX +1 — 3NX, 3NX +1 — n. This numbering system 
Is chcsen so that the resulting variables in the computer listings 
are identified through double subscripts only. The first subscript 
Identifies the line along which the variables are calculated vhilt 
the second subscript indicates the position along the line. In this 

work the lines are constantly spaced with spaclngs h . h and h . 

X y z 

This is done purely for convenience. The advantage of uniform line 
spacing is that the resulting ordinary differential eqviations can be 
solved more easily than those that are derived fron non>unifonn line 
spacing. 

The equations of 3>D elasticity are the three couplad particl 

differential equations (2.5-2. 7). In this case solutions for the 

dependent variables are possible only at points where the three sets 

of parallel lines intersect. These points are usually called nodes in 

a discretized body. This limitatic is the result of coupling among 

the equations, which makes the particular selection of ordinary 

differential equations valid only at the nodes. The equations are 

developed for a crack specimer as shown in Figure A. Sin«.e we have 

a three fold symmetry, only one quarter of jhe specimen is considered. 

One should also note the numbering of the various planes. 

2.2.1 Reduction of the First Navler-Cauchy Equation and Associated 
Boundary Conditions for the Cracked Specimen 

For the solution of equation (2.3). the lines parallel to x axis 
in Figure 3 are considered. The x directional dlsplaccmantu of points 
along these lines will be denoted as \x^ u^, — u^^. We define vj^^, 
v| 2 . VI 3 . . .vj^as the derivatives of the y directional displace- 
ments of the same points on these lines with respect to y and v|^, 



15 


ORIGINAL PACK IS 
OF POOR QUAL'TY 


v] » w| - • • • w| p as the derivatives of the z directional displace- 
ments of the same points on these lines with respect to z. These 
displacements and derivatives can then be regarded as a function of 
X only since they are variables upon lines which are parallel to the 
X axis. 

By the above given definiton, the ordinary differential equation 
along a generic line ij in figure 5 may be written as, 


d^u 


IJ 


dx 


2 2 

+ c (i-y + 

1 


d * I ^ d • I 

"2 


(2.15) 


Introducing 5- point finite differences, the partial derivatives of 
u with respect to y and z along the x directional line (ij) of 
figure 5 can be written as f.^xlows. 


' 2^i1 2 

ay ^ 12h^ 




30uij + 


(2.16) 


and 


(~) 


12h 


r ■ “1.3+2’ 


(2.17) 


On substituting equations (2.16) and (2.17) into equation (2.15), the 

general equation along Interior lines is obtained. Thus, ... 

(.2 . lo) 

^ ^^^1^ ^"V2.j'^ ^^Vi,j ■■ '’1+2, j^ 


dx 


<-“1,3-2 <^“1.3-1 - ’““13 * ““1.3+1 - “1.3+2” "=2 dJ <“ + "’13 
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The rearrangement of Che above equation leads to, 
C. 4C 

1 . f i . 

T~ 4 ■*“ <-T- + 

2h 


d^u 


1 -1 ^^1 . 

T “i-2 1 r Vi i “T^ 

dx" 12h‘^ ^ 3h^ ^ 2h^ 2h 

y y y s 




4C. 


3h' 


^1 ^1 

y act 


s 


(2.19) 


Similar differential equations are obtained for Che other displace- 
ments of the points on the x directional lines. Since each equation 
has the terms of the displacements of the points on the surrounding 
lines, these equations constitute a system of ordinary differential 
equations for the displacements Uj, . 

The equation (2.19) is applicable only to interior lines. For 
boundary surface lines and lines adjacent to the boundary surface 
lines, the difference expressions for the second derivative will in- 
volve imaginary lines outside the boundary. As shown in figure 6, 

f y f V 

the imaginary lines below the plane V are designated as 1 , NY+1 ' 

. . . (? - NY+1)^^. The imaginary lines lying on the planes adjacent 
to plane VI are numbered as 1^^* . . . and 1^^* . . .NY * respectively . 

Similar numbering procudure is followed for the imaginary lines lying 
on the planes, which are located adjacent Co planes II and III re- 
spectively. The arrangement of the planes V, VI, II and 111 is shown 
in figure 4. Since three-dimensional elasticity problems have three 
boundary conditions at every point of the boundary surface and a 
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second order ordinary differential equation needs only two conditions* 
the shear stresses at the boundaries are used to eliminate the 
imaginary lines outside the surface. While the condition of the pre- 
scribed normal traction or displacement will be enforced through the 
constants of the homogeneous solutions. The method used to derive 
the equations for boundary surface lines and line adjacent to it, is 
described in detail in the Appendix A. Appropriate difference equa- 
tions are used to express the partial derivatives. For example » to 
write the difference equations for the partial derivatives with re- 
spect to z for the lines t - NY+1 through 2, a backward finite dif- 
ference formula is used due to the lack of sufficient grid points to 
write central difference equations. On following the procedure des- 
cribed in Appendix A, the differential equation for line 1 can be 
written as. 


^ [-20u^ + 17 m2 + '>U 3 - 


dx 12h 


+ ^ [-30i. + 30u - 2u^ 

12h‘ 


NY+1 


^ lie, , 

_ p d , • • X 1 dv I 

^2 dx ^1 * 6h dx ‘l 


U.20) 


similarly the differential equation for the line 2 - 2NY+1 is 


^ - 2NY+ 1 ^1 , _ , ,, 

2 - 2NY+1 ^“2 - 2Nl['+2 


dx 


12h 


^"2 - 2NY+3 ■ “2 - 2N’Y+4) + — ^ X 

12h" 


z 
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- 4NY+1 ■*' - 3NY41 " - 2NY>1 ^ - NY-1^ 


'2 dx 


(v + w) 


8, - 2NY+1 


^1 

6h dxk - 2NY4-1 

y 


6h dx ' il - NY+1 

z 


( 2 . 21 ) 

Finally, Introducing matrix notation all the ordinary differential 
equations along the x directional lines are expressed in the form. 


- [K,] (u) + ^ {r(x)} 
dx 


fxl 8.x88.xl £.xl 


( 2 . 22 ) 


where the matrix [K^] and the column vectors {u} and (r(x)} are 
given below. 


(Kijl 

NY X NY 

"l2> 

tn K NY 

1N„1 0 

NY X NY 

0 



0 

0 

^ 1 

iKjjl ! IHjjl 
NY X NY NY' X NY 

1 

i«ji' i 's*' i “ 

MY « NY . NY « NY 

0 

— 




"al' 1 

|N„1 

|K„) 

1*^2^ j 

0 1 

0 1 

NY X NY 

NY X NY ' 

NY' X NY ^ 

NY X NY ' 

1 

NY X NY 


i 

0 





s. 

e 

0 

1 

° i 


Uj,l 

|K„1 

"‘J.' , 

IK„I 


1 

1 

NY X NY 

NY’ X NY 

NY X NY 

NY X NY;NY X NY 

0 

0 

0 


tK„l 

|K„I 

INj,) 




NY X NY 

NY X NY 

NY X NY 

NY . NY 

0 

0 

0 

1^2,1 


^*NZ,NZ-l^ 

*SiZ,NZ^ 

_ 



NY X NY 

NY X NY 

NY X NY 

NY’ X NY’ 

_ 


(2.23) 
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[K^2l 


- I V [I] 

z 




[ 1 ] 


6h‘ 


[K2i1' 


3h 


[I] 




12h 


[I] 


%Z,NZ-2' “ " ^ 
Jn 


II] 


^^Z,NZ-1^ 


17C, 


12h 


(I) 


Matrix [1] is an indentity matrix of order IIY x NY 


Note that is 


1 > 2v 
2^1 - v) 
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{u} 

Jlxl 



NYxl 



(x)} 



(x)} 


(r(x)} = 

8,xl 


{u> 


NZ-1 





where the partitioned column vectors are 


" 

“l 


'^NY+l 

L 

‘2 


Vy+2 

1 

1 

1 

1 

1 

U 

NY~1 

NYxi 

i 

1 

1 

1 

1 

1 

1 

1 

“2NY-1 

“ny 


'^2NY 


(2.2A) 


i 

\ 

i 
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(r^Cx) } 
NYxl 



. . llC.v 

"'ny+i - sr - 


y NY+1 


. . C.v 

(v 4- v) 4- — ^ I 

''-'nY+ 2 ^ 6h I 

y NY41 


‘^2 ‘^**)nY +3 


S «■' + ''>2NY-2 


(v 4> w) - — iL. I 
NY-1 6h I 

y 2NY 


. . llC,v 

(v 4- w) '4 

^2NY 6h 


y 2NY 
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llC-v 

I I 


‘'z ''^l-2NY+l - -6h 


y 4-2NY+1 ^*^2 i-J 


C V C.w 

‘=2 '''*“>i-2OT+2 ■'jrl -sri 

y a-2l^Y+l ®"2 H-NY+2 


S ''^i-2NY+3 ~ 6h 


I i-NY+3 


<^2 “>!- OT -2 - 6 K 


2 1-2 


. . C V Cw 

S<'' * ">2-NV-l - 


llC,v C w 
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2.2.2 Reduction of the Second Navier-Cauchy Equation and Aasociated 
Boundary Conditions for the Cracked Specimen 

For the solution of equation 2.6, the lines parallel to y- axis 

in figure 3 are considered. The y directional displacements of points 

along these lines will be denoted as V 2 , v^, .... v^. We define 

uL,uL, uL u as the derivatives of the x directional dis- 

1 2 3 m 

placements of the same points on these lines with respect to x and 
w|j^, wl^.wl^t • . derivatives of the z directional dis- 

placements of the same poj^^s on these lines with respect to z. These 
displacements and derivatives can tu^n be regav'^ed as functions of y 
only since they are variables upon lines which are parallel to the y 
axis. 

By the above given definition, the ordinary difterential equa- 
tion along a generic line ij may be written as. 


dy dx dz 






(2.25) 


Introducing 5- point finite differences, the partial derivatives of 
V along the y directional line iJ can be written as follows. 




12h 


^ "'VlO - ’<^13 ^ ‘‘V 1.3 - ''1+2,3^ 

(2.26) 


and 


3"v 


12h* 


<- ''1,3-2 “''1.3-1 ■ ’“^13 “’'i.3+1 ■ ''1.3+2 


(2.27) 
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On substituting equations (2.26) and (2.27) into equation (2.25) the 
general equation along interior lines is obtained. Thus, 


d^v 


dy 


^ + C. [ — ^ (-V -j . + 16v - 30v + 16v . - V ) 

2 1 j^2h^ il 1+1.3 i+2,j 


12h 


T~ 4 0 + l^'^- 4 1 “ 30v + 16v, -V, . , -)] 

2 i,j-2 i,j-l ij i,3+l i,j+2'^ 


^ ">ij 


(2.28) 


The rearrangement of the above equation leads to 


dv.. C, 4C, 5C, 5C- 

ij 1 1 j_/lj^l\ 

dy^ 12h^ 3h^ ^ 2h^ 2h^ 

X X X z 


4C, 


AC, 


2 i,j+l i.j+2 3^2h^ i»3"2 3^,2 i,j-l 


3h 


• • 


+ C j- (u + w) 

2 dy ij 


- V + ■ - ■»— V 

3^2 ^+1.3 


12h 


2 i+2,j 


(2.29) 


Similar differential equations are obtained for the other displace- 
ments of the points on the y directional lines. Since each equation 
has the terms of the displacemnts of the points on the surrounding 
lines, these equations constitute a system of ordinary differential 
equations for the displacements , V2 • • • v^. 
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The equation (2.29) is applicable only to interior lines. For 
the boundary surface lines and lines adjacent to the boundary surface 
lines » the finite difference expressions for the second derivative 
will involve the displacements on Imaginary lines lying outside the 
boundary as shown in figure 7. The shear stress at the boundaries is 
used to eliminate the displacements on these Imaginary lines near the 
surface. While the condition of prescribed normal traction or dis- 
placement wil be enforced through the constants of the homogeneous 
solutions. 

A detailed description of the method used in deriving the equa- 
tions for the lines near the boundary, can be found in the Appendix A. 
Tlie use of fictitious lines and their subsequent elimination using 
the conditions of symmetry and prescribed shear stress at the surface 
leads to the following equations for line 1, 


^ + 17Vj^2+1 ‘**2NZ+1 " ''3NZ+; ^ 

dy 12h 

X 


+ s- [-30v 4 32v - 2v ] 

12h 1 / J 


^ lie, 

S d-y ^ ~ 'l 

6h 

X 


(2.30) 


on following the same procedure the ordinary differential equation 


for line m can be %rritten as. 
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— ^ + ^ [ “V + Av + 17v - 20v ] 

, 2 10.2 m-3NZ m-2NZ m-NZ m‘ 

dy 12h 


+ jp- [-V - + 4v ^ + 17v . - 20v ] 

-_.2 m-3 m-2 m-1 m 

12n 


, lie. , lie, , 

„ d-*v . Idui . Idwi 
C_ — (u + w) + 7T — T~ 7T — T~ 

2 dy m 6h dy m 6h dy m 


(2.31) 


The comparison of equations (2.30) and (2.31) with equation (2.28) 
shows the presence of additional terms on the right-hand side. These 
terms appear here due to the consideration of ioiaginary lines outside 
the boundary surface. 

Finally, introducing matrix notation, all the ordinary differen- 
tial equations along the y- directional lines are expressed in the form 


^ • 'S' (''1 * fy 

dy ^ 

mxl mxm mxi mxl 


(2.32) 


where the matrix the column vectors {v} and {s(y)} are given 


below. 
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NZxNZ 








— 

5Cj 31Cj 

“T ~~T 

2h , 12h'^ 

z X 

2 

s 

■ ~r 

6h 

2 

0 

0 

0 

0 


31Cj 31Cj 

I2h^ 12h^ 

2 ) 



0 

0 

0 

z 


12h^ 

2 




s 



‘S 

s 

0 

0 

12h^ 

2 

2 

2h* 122i 

2 X 

2 

12h^ 

2 



0 

K 

% 

*k. 

V 

V. 

V 

r . 

% 

^ 

V 

*■ V 

V 

0 

0 

0 








12h^ 

£ 



2 

12h^ 

2 

0 

0 

0 

12h^ 

2 

‘S 

31Cj^ 31Cj 






. s 


5C^ 31C, 

- 0 ■ • 

• -0 

0 

12h^ 

2 

3h^ 

2 

12h^ 

2 

3h' 12h^ 

2 X 
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iKi^l 


IK13I 

[K21] 


17C, 


m 


12h 


3h' 


[I] 


m 


12h 


4C, 


3h‘ 


m 


Matrix [I] is an identity matrix of order NZ x NZ. Note that 


Is 


1 - 2v 
2(1 - V) 


“ 


- 

(v)^ 


{S3(y)} 

(V), 


{s2(y) ) 

1 

1 

1 

(s(y) } 


t 

1 

mxl 


1 





(s^(y)) 

■ m 


- 


(2.34) 
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where the partitioned column vectors are 


- 


“ « 



'^NZ+l 

^2 

1 

1 

1 

1 

1 

‘v>2- 

'^NZ+2 

1 

1 

1 

1 

1 

'^NZ-l 

NZxl 

1 

'^2NZ-1 

'"nz 


'^2NZ 

- 


- 



w _ 


r 


'"in-2NZ+l 


'^m-NZ+l 


''m-2NZ+2 

1 


'^m-NZ+Z 

{v} = 

1 

1 

1 


1 

» 

t 

1 

NZxl 

V 

m-NZ-1 

NZxl 

'^m-l 


''m-NZ 


V 

m 


«• ^ 
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(sj^(y)} 

NZxl 


. . lie u 

S < “ + “>1 - -6h- 


X 1 


• » 


lie u 

(u + w >2 - -gjp- 


X 2 


llC^u 

^2 ''^NZ-2 ■ “6h“ 


X NZ-2 


• • 


llCjU C w 

C, (u + W) r^l - 

^ NZ-1 NZ 


. . iic» 

C,(u + w) --^ + -7T^ 

2 NZ 6h^ 6h^ 
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r 


C, (u + 6h 


• • 


X 1 


{s2 (y) } 

NZxl 


• • 


CiU 

^2 ”^NZ+2 ■"ST', 

X 2 


C u 

• • . i ] 

C- (u + w) „ + r^l 
2 2NZ-2 6h^ 


• • 


C^u C^w 

‘^2 "’2NZ-1 ’ ^'zNZ 


« • 


C^u llC^w 

S <''^“’2NZ^6Tj,/-6h7l,,, 
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CiU 


NZxl 


♦ • 


V, 


S “’m-aNZ+Z 6h^'„.Nz^.2 


S“, 


V, 


^2 ^^m-N2-l ■ 6h ' 6h 

X m-l z 


• • 


C.u llC.w 
^2 ■*■ ''^m-NZ “ “6h“l 

X m z 


m-NZ 


m-NZ 
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2.2.3 Reduction of the Third Navier-Cauchy Equation and Associated 
Boundary Conditions for the Cracked Specimen 

For the solution of equation (2.7), the lines parallel to z- axis 

in figure 3 are considered. The z- directional displacements of points 

along these lines will be denoted as w, , w.. , w. ,...., w . We de- 

1 2 3 n 

• » *1 *1 • 

fine uL,uL, ul*, ...,u as the derivatives of the x- directional 
'1 2 3 n 

displacements of the same points on these lines with respect to x 

and • • •» as the derivatives of the y- directional 

displacements of the same points on these lines with respect to y. 

These displacements and derivatives can then be regarded as functions 
of z only since they are variables upon lines which are parallel to 
the z axis. 

Following the procedure as before, the ordinary differential 
equation along a generic line ij may be written as, 


,2 
d w, . 


— ^ ^ (-W. - . + 16 w. . . - 30 w, . + 16 w . 

^ 12h 2 ^-2-3 ^-1-3 ^3 "+1-3 


- w 


i+2,j) 


+ ;; (-W. . + 16 w. - 30 w + 16 w - w . . ] 

^2^ 2 1,3-2 i,j-l ij i.J+1 i,j+2 


■ '^2 d 1 + ''>11 


(2.35) 


Similar differential equations are obtained for the other dis- 
placements of the points on the z- directional lines. Since each 
equation has the terms of the displacements of the points on the 
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surrounding lines, a system of ordinary differential equations is ob- 
tained for lines w, , w* . . w . The equation (2.35) is applicable 
1 Z n 

only to interior lines. For boundary surface lines and lines adjacent 
to the boundary surface lines, the finite difference expressions for 
the second derivative will need the use of the fictitious lines. The 
displacements on these fictitious lines is eliminated with the help of 
shearing stresses on the free surfaces. 

As shown in figure 8, the z- lines 1 through NX are split into 
two sets. The first set consists of lines 1 through NXC lying on the 
plane Vtb). For these lines a fictitious line is used to write the 
finite difference equation, which is subsequently eliminated by using 
the shear stress boundary condition on this face. The second set con- 
sists of lines NXC + 1 through NX lying on the plane V(a). The w dis- 
placements are symmetric with respect to this plane. Therefore, 
symmetry consideration helps in eliminating displacements on the fic- 
titious lines. A detailed description of the method used in deriving 
the equations for the lines lying on the boundary and adjacent to the 
boundary, can be found in the Appendix A. For example, the ordinary 
differential equation for line 1 is 

d^w. C 

2 T '^l ^ ”nx+ 1 ^ ''ZNX+l “ ''SNX+l^ 

dz 12h 

y 

Cl 

+ s- [-20 w^ + 17 Wj + A w- - w, ] 

12h 


X 
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^ lie, . lie, , 

^ d 1 du I 1 dV| 

^2 dz '"H ■ 6h dz 'l “ 6h dz'l 

X y 


(2.36) 


For line NX - 1 the equation can be witten as; 


+ — r-30w +32w -2w 1 

2 ,0^2 ' NX-1 2NX-1 3NX-1^ 

dz 12n 


~,2 ^"''nX-3 ''nX-2 “ ''nX-I ''nx^ 

Izn 

X 


„ d , . , . 1 du I 

^2 dz '"^NX-1 ~ 6h dz' 

X NX 


(2.37) 


Finally introducing matrix notation, all the ordinary differential 
equations along the z-directional lines are expressed in the form. 


2 

= [K ] {w} + {t(z)> (2.38) 

,2 z dz 

dz 

nxl nxn nxl nxl 

where the matrix [K ] and the column vectors {w} and {t(z)} are given 

z 


below 
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rx,i 

nxn 


"'ll! 
NX X NX 

iKjj] 

NX X NX 

IK13) 

NX X NX 

(K141 

NX X NX 

0 

0 

0 

[K21I 1^22! 

NX X NX NX X NX 

^ 1 

(Kjil 1 [Kj^] j 

NX X NX 1 NX X NX 

0 

0 

0 

[ , 

1X2^1 [K2jI ; [K2jI 

1 1 

' I i 

NXxNX NXxNX NXxNX^NXxNX 

(K241 

NX X NX 

i 

0 1 0 

1 

i ^ ^ ^ i * • ^ 

0 ^ 1 ^ 

1 ^ ^ r V 1 ^ • * 

' 1 I' 


0 

0 

0 

1 

IK24I 

NX X NX 

lK,j) 

NX X NX 

'■'33' 

NX X NX 

tK^il I [KjJ 

1 

NX X NX'NX X NX 

0 

0 

0 

IK2^1 

NX X NX 

IK23I 
NX X NX 

IK22) 

NX X MC 

f'21> 

NX X NX 

0 

0 

0 

f'241 

NX X m 

Sy,NY-2^ 
NX X NX 

^Sy,ny-J 

NX X NX 

^Sy.ny^ 

NX xNX 


( 2 . 39 ) 
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where 


IK 


12 


(K^^l j 101 

-----T “1 

[01 :[K^2l 


17Ci 

'^^NXC.NXC 

y 


8C, 


IK12I 


[I] 


3h 


NX-NXC, NX-NXC 


[K,3l - 


[K^jl* j [0] 

[0] |Fk^ 31" 


* 1 

1*^13! ■ ■ ^^Inxc.nxc 

y 


[K13I 


.. 


[Ki,] - 

NXxNX 


6h^ 

y 


' 

: 10] 

• 

[0] 

1 [0] 


[K^,] 


T ^^Inxc.nxc 


12h 


4C^ 

^‘^2il " ’ TTF ^^Inxxnx 

3h 

y 


“ ,,, 2 ^^‘NXxNX 
12h 


y 
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^‘SlY.NY-2^ “ " rr f^^NXxNX 

jn 

y 

^^.NY-1^ " " ^^WnX 

y 
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lie u C V 

4 ^>1 - "e^T"' 6 t1 

X NX+1 y 1 


C u C V 

c_ (u + v)^ + i + rr~l 
2 2 ^\^^x+l 


C^v 


C (u + V) + 


y 3 


C2 (u + v)nxc + 

y NXC 


{t^Cz)} • 
lOCxl 


■^2 <“ ^ *>NXC+1 


S * ^^NX-2 


CiU 

2 ^ 2NX 


C, (u + v) + 
NX 


lie u 

1 

6h 


2NX 


Nxexi 


(NX-NXC)xl 
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NXxl 


llCjU 

'2 2“'+! ■ “^^KX+l 


C^u 


C- (u ^ V>2J^ 2 ■*■ gfi“' 

2 2NX+2 2NX+1 


^2 ■*■ ^^2NX+3 


C2 (u + 2 


^2 '^W-l ~ 6 h ^ 


C,u 
-! 

X 3NX 


lie u 


X 3NX 
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^2 n-2NX+l 


llC^u C^v 

6h ' I n-2NX+l ^^n-NX+1 
X y 


CjU C^v 

^2 '’■ ''^n-2NX+2 '‘' 

X n-2NX+l y n-NX+2 


C 2 (u + _ 1 j 

, ^ y n-NX+3 


C2 (u + I _2 


S + ''>„-,x-l - - 6h^ 


llC^u C^v ^ 

^^2 "’n-NX ■*■ 'uT'n 

X n-NX y 
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NXxl 


C 2 (0 + v)^ 


S ^>n-NX+2 


llC,u 

^ 1 1 

llC^v 

6h ' ' 

6h 

X n-NX+l 

y 

C,u 

11c, V 

-i-l + 

1 

n-NX+l 

6h 

y 


llC.v 


=2 <“ * *>n-NX+3 "Sir 


y n-NX+3 


llC,v 




C, (u + v) + 
dC n 


6h 

i 

y 

n-2 

C u 

llC,v 

_ 1 I 

+ ^ 1 

6h ' 

6h I 

X n 

y n- 

llC.u 

llC,v 

1 1 

+ ^ 1 

6h 1 

6h 1 

X n 

y n 


n-NX+l 


n-NX+2 
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2.3 Solution of Differential Equations with Constant Coefficients 
The set of I second order differential equations represented by 
equation (2.22) can be reduced to a set of 2? first order differential 
equations by treating the derivatives of the u's as an additional set 


u 


Jl+1 



“ s ,+2 



etc. 


The resulting 21 equations can now be written as a single first 
order matrix differential equation 


dx 




dR(x) 

dx 


(2.41) 


where U and R are column matrices of 2£ each and is 21 x 2H matrix 
of constant ..oef f Icients and are written as follows. 


U - 

2HX1 


Hxl 


dU]^ 

dx 

du2 

dx 


‘‘"i 

dx 


Axl 


(2.42) 
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R - 

29x1 



{r(x> i 


£xlj 


(2.43) 


and 


”1 

21x21 



.L « 




(2.44) 


In the similar manner the second order ordinary differential equation 
given by eqution (2.32) can be transformed to a first order ordinary 
differential equation and it is given by, 

^ - A V + (2.45) 

dy 2 dy 


where matrix A 2 and column vectors V and S(y) are 




^1 




^2 

1 

1 

1 





mxl 

V - 

• 

M m 


2mxl 


dx 

dx 

dv 

m 

dx 

m « 

mxi 

. 


(2.46) 
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S(y) 

2mxl 
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and 


10] ' U] 

mxm » mxm 


imxzm 


[K ] ' [0] 

y mxm ' mxm 


(2.47) 


(2.48) 


For the z- directional lines the second order differential equa- 
tion given by equation (2.38) is also transformed to a first order 
ordinary differential equation. The transformed equation is. 


dw’ 

dz 


W + 


dT(z) 

dz 


(2.49) 


where matrix 


and column vectors W and T(z) are given below. 


w’ 

2nxl 


- 


• 


'^1 

”2 

1 

1 



w 

nxl 


L ^ J 



dWj^ 

dx 

dw^ 

dx 

• 

• 

dw 



n 

dx 

nxl 


• 




* 


(2.50) 
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( 2 . 52 ) 
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2.3.1 Solution by Recurrence Relations 

The system of ordinary differential equations (2.41), (2.45) and 
(2.49) can be solved by any number of standard techniques. An algorithm 
is derived for obtaining the solution by using a recurrence relation to sweep 
through, from one boundary to the opposite one, once, to obtain the 
missing boundary conditions and then to sweep through a second time to 
obtain the complete solution everywhere. These equations are not in- 
dependent of each other. They are coupled to the other equation through 
the vectors R, S, and T. Therefore, the solution of the equations can 
be achieved by an iterative technique in which the solution for two 
sets of different directional lines is used to evaluate the coupling 
terms for the third set of directional lines. Then the complete solu- 
tion for this set of directional lines is achieved through the recurrence 
relation method described below. 

Consider the matrix equation 


r = A'' + 

dx dx 


(2.53) 


Let the x lines be divided into n intervals (not necessarily equally 
spaced) with h^ the spacing between the (i-l)£l^ and i th nodal points. 


Then 


U. 

1 



^i^i ^i-l^i-1 ^i " ^i-1 

n • 1 


(I - ^ Apu^ 


(I A y A (R^ - R^_p 


(2.54) 


Let 


64 


ORlGf^■Al 

OF POOR QUALH Y 


or 


-1 ^1 
. (I . y A^) (I 4 ^ 


+ (I - Y A^)"^ (R^ - R^_j) 


(2.55) 


or 


where 


”i ■ h h-i * ”i 


-1 

(I - Ai) ^ (I + y A^_^) 


- (I - y A^)‘^ (R^ - R^_^) 


(2.56) 


U, - 4 


(2.57) 


Then 


or 


®i “i *h'h «i-i + «i 


■ h ^®i-i"i * * «i 




Since Uj^ Is abritraty 


“i - h “i-i 


F. “ L F + M 
1 1 1-1 1 


(2.58) 


65 


ORIGINAL PAGJI 
OF POOR QUAUflf 


For station 2, we can write from (2.57) 




and from (2.55) 




^2 " ^2 “ ”2 


(2.59) 


All the and F^ can now be determined from 2.58, starting with (2.59) 


At the last station, i - n, we have 


U = D U, + F (2.60) 

n n 1 n 

Now parts of the vectors and are known (usually half the boundary 

conditions will be given at each end). Let us write these vectors as 


• 

“l,k 


f 

c 

“l.U 

II 

c 

U 

n,u 

■■ W 




(2.61) 


where k represents the known part of the vector and u the unknown part. 

Let us further partition D and F as 

n n 


r 1 


- "> 

^ «2 


Bn 


F 

1 


. n 


do 01 , 



3 4 


^2 



. 
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So that we can write 


• - 


r 

•i 


► ^ 


• • 

u 1 . 

n,k 

■ 


^2 


"l.k 

+ 

®1 

$2 

U 

n,u 


“3 



"l.u 


m m 


— 

• 


■» m 


ft a 


Then 


n,k 


a, U, . + a, U. + 6. 
1 l,k 2 1 ,u 1 


(2.63) 


or 


Ui - a~^ 

1 ,u 2 


^^n,k ■ '^1 ^l.k " ^1^ 


(2.6A) 


The vector Is now completely known and the complete solution obtained. 

It is apparent that if the first part of the vector is the unknown 

part and/or the first part of is unknown, then the appropriate relation, 

taken from equation (2.63), must be used to determine U. 

i ,u 

Thus, suppose the first half of is unknown and the second half 
known, so that 


‘'1 ■ 


r n 

“l.u 

I'l.k 


Then from eq. (2.63) 


U . - a, U- + a_ U. , + B, 
n,k 1 l,u 2 l,k 1 


and 


U, - a~^ (U . ” a» U- - BJ 
l,u 1 n,k 2 l,k 1 


(2.65) 
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similarly if tht- first part of is unknown and the second part is 


known. 


then 


U 

n 


“ 

U 

n.u 



j 


and 


n,k 


= a. 




a, U, + 
A 1 ,u 


and 


U 


1 


»u 





U 


l,k 


Bj) 


and finally if 


_ ■ 

u. 


U 

l,u 

and U *= 

n.u 

U , 

n 

U , 

l.k 


n.k 

M * 


• 


n.k 


a. U, + a, U. . + 

3 l,u A l,k 2 


and 


U, = a'^ (U V - V - 

l,u 3 n.k A l,k 2 


( 2 . 66 ) 


(2.67) 


To summarize: 

1. The vectors and are known (R obtained by iteration). If 
the spacing between lines is maintained constant, then E A is a con- 
stant matrix. 


2. Calculate the and from (2.56). 
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4. Partition the laatrlces and as indicated by equations 
(2.62) and (2.63). Note that if there are M lines then the vector U 
will have 2M elemenva and there will be 2M boundary conditions, divided 
between stations 1 and n. This will determine the length of the known 
and unknown vectors and thus the dimensions of the a*s and B*s. 

5. Calculate the unknown vector U. from (2.64) or (2.65-2.67). 

i. ,u 

The complete vector is now known. 

6. The displacements at any station are now computed from (2.57). 

As noted above if the spacing between lines is constant, l.e., h , h , 

X y 

and h are constants, then A is a constant matrix. It then follows 
z 

from (2.56) that - L is a constant matrix (if h^^ is ilso constant). 
Then equations (2.58) read 


LD 


i-1 


LF 


1-1 


+ M. 


Note that although varies along the lines, it is only a column matrix. 

2,3.2 Incorporation of Prescribed Boundary Conditions of Normal 

Stresses or Applied Displacements in tha Recurrence Relations. 

The second order differential equations given by equation (2.22) 

for the X-lines involve a two point boundary value problem. The first 

half of the boundary conditions are given on face IV and the remaining 

boundary conditions are given on face I as shown in Figure 4. In the 

case of an edge crack specimen they are 


0 


XX 'IV 


0 


( 2 . 68 ) 
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and o L ■ 0 (2.69' 

The subscripts IV and I refer to the boundary planes as shown in 
Figure 4. With the help of stress-displacement relations the equations 
(2.68) and (2.69) are reduced to the following equations respectively. 


du I 
^'IV 


V 

1-v 


(v + w) 


IV 


(2.7C) 


and 


du I 
dx I 


V 

1-v 


(v + w) J 


(2.71) 


Since the derivatives of v and w will be known on all the points on 
the face IV and I, a vector of {U^ k^£xl ^^'n k^lixl formed by 

assembling the derivatives of u in the proper order, as designated in 
Figure 6. These vectors will constitute the lower half of the complete 
vector used in the recurrence relation method. The complete set of 
vectors are 


and 


= 

2{.xl 


U = 
n 

2ixH 


(U 

{U 

{U 

{U 


1 ,u^£xl 
1-k^Jlxl 

n,u^P,xl 
n,k^ .'ixl 


r-.72) 


(2.73) 


As described in the previous section using equation (2.66) the remaining 
half of the vector can be calculated and a forward sweep is made to 
obtain the solution at every grid point. This procedure presents a 
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problem because bt the boundary conditions are Neumann's type and 
due to this the solution can be determined only upto an arbitrary con- 
stant. The presence of rigid body displacements makes matrix cx^ in 
equation (2.66), a singular matrix. To suppress the rigid body dis- 
placement, the displacement Is set to zero. This leads to a 

rearrangement of matrices and in equation (2.66). Once 
{U. ^ is determined by using the rearranged equations, a forward 

sweep is carried out to determine the solution at all the grid points. 
Th* known value of saved before the sweep is carried out. 

To maintain the consistency the newly calculated value of 
replaced by its old value. 

The second order differential equation given by equation (2.32) 
for the set of y-dlrectional lines is again a two point boundary value 
problem. Half of the boundary conditions are prescribed on face V 
and the. remaining half of boundary conditions are given on face II, 

The face V and face II are the boundary pxanes as shown in Figure 4. 

It also shows the cracked face designated by the plane V(b) , This 
leads to two kinds of boundary conditions for the points lying on the 
plane V. The face V(b) is e. traction free surface, therefore all the 
y- lines, starting from this face will satisfy the following boundary 
condition. 


^yy'v(b) " ° 


(2.74) 


Due to the symmetry of v-displacements with respect to the plane V(a), 
all the lines starting from this face will satisfy the condition. 
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= 0 


(2.75) 


using the stress-displacement relations equation (2.74) can be reduced 
to 


dv I _ V ^ du ^ 

V(b) ' ' V(b) 


(2.76) 


The equations (2.75) and (2.76) are assembled into a vector form by 
following the sequence of numering of the y-directional lines as shown 

T II 

in Figure 7. The assembled vectors are and Using these 

vectors, the vector at face V is written as 


"i 

2mxl 


l,u 


l.k 


l,k 


^l,u 

L J 


where v' and v" represent the unknown parts of the vector V. . 

Rest of the boundary conditions are obtained by utilizing normal 
traction at the face II. As shown in Figure 2, in the case of tensile 
loading we have an applied normal traction 0 on this face. Using 
the stress-displacement relation, we can write. 
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dv I ^ r_a_ (1 + v) (1 - 2v) 1 _ V , du 

dy'^^ E (1 - v) (1 - V) Mx jj 

(2.77) 

The equation (2.77) is arranged to form the vector V , . The vector 

n, A 

V at face II can now be written as 
n 

(2.78) 

V represents the unknown par*" of the vector v‘ Following the 
H ^ u 

recurrence relation method, an equation similar .o the equation (2.57) 
is formulated and given by 

V^-DnVi-fF^ (2.79) 

The evaluation matrices D and F is described in the previous sec- 

n n 

tion 2.3.1. The determination of unknown vectors v’ and V*' differs 

l,u l,u 

slightly from the procedure given in the section 2.3.1. The details 
of the procedure used to calculate vector can be found in Appendix 
B. Once the complete vector is known, using the recurrence relation 
method, the displacements v and v are obtained at all the grid points. 

Tht. second order ordinary differential equation for the set of 
z- directional lines given by the equation (2.38) is a two point boun- 
Gary value problem. As shown in Figure 4, half of the boundary conditions 
are obtained by setting ill the w- displacements equal to zero on the 
plane VI due to symmetry. This leads to 
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“ 0 (2,80) 

nxl 

on the boundary plane III» the prescribed normal tractions are zero. 
This leads to a set of equations for each grid point on this plane. 
These equations In vector form can be written as 


W - 

n.k ^dz^ 

nxl 


nxl 




(2.81) 


once again the vectors W. and W are formed as 

1 n 


and 




2nxl 


W 


l.k 


W 


l»u 


W 


w 


n,u 


W 


n.V. 


(2.82) 


the general equation connecting the vectors and is 


W « D W- + ? 
n n 1 n 


( 2 . 


The matrices and can be partitioned as shown in equation (2.63). 
Appropriate matrix equations are extracted to solve for The 

complete solution for w- displacements and its derivative w is obtained 



74 


ORIGINAL P/ 

OF POOR QoAl.TY 


for all the grid points by carrying out a forward sweep operation, 
described In section 2.3.1. 

2.4 Loading Idealization for Compact Tension Specimen 

One of the major objectives of this Investigation is to study the 
compact tension specimen. Figure 1(a) shows a standard compact tension 
specimen and the loading on It. The Idealized model of this specimen 
used In the present study is shown In Figure 1(b). The pin loading 
applied on the apec:<men Is approximated by a parabollcally varying 
shear type loading. The resultant force due to the applied shear load- 
ing Is maintained equal to the actual applied pin load (25). 

This type of prescribed shear loading can not be Incorporated into 
the solution of the ordinary differential equations in which the noirmal 
traction boundary conditions have been used. The equations for y- lines 
are slightly modified to Include this type of traction boundary condition. 
The y- directional lines whose equations are directly affected by this 
are 1, 2, . . NZ and KZ + 1, NZ + 2 . . 2NZ. 

The traction a at the boundary face IV is 
ny 


a 


ny 


IV 


-T(y) 


where t(y) is the prescribed surface . 
direction is expressed as 

T(y) - 4 T^(^) (1 - Z) 

T is the maximum value of prescribed 


(2.84) 

traction. Its variation xn y- 


(2.85) 


surface traction at y 


h 

2 


and h 


is the semi length of the specimen as shown in Figure 1(b). The 
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inclusion of this leads to a new equation for the y- directional line 
1. The modified equation is 


d^v. 


dy 


12h 


[- 20v^ + 17vsz+l + 


— V 1 

3NZ+1^ 


12h 


^ [- 30v^ + 32V ^ - 2V^1 


„ d 1 di 

2 dy 1 oh dy 


I [Illl] 

^ G ' 


( 2 . 86 ) 


where G is the shear modulus of the material. Similarly the equations 

for lines 2 through NZ can be written. The equation (2.86) is the 

llCi (y) 

same as the equation (2. 30) with the additional term of - [ q 1 » 

which appears due to a non-zero prescribed surface traction. The 
modified equations for line 2 through NZ also include this additional 
term. The other terms for these equations remain the same as before. 


The ordinary differential equations for the line NZ+1 through 2NZ are 

^1 x(y) 

modified by adding additional term of g — 


to the original equations. 


The details of the method used to derive these equations can be found 
in the Appendix C. 


Introducing the matrix notation, all the ordinary differential 
equations along the y- directional lines are expressed in the form. 
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Y 


{v} + ^{s(y)} + (8*(y)} (2.87) 

mxm mxl mxl mxl 

The equation (2.87) is similar to the equation (2.32) with an addition 
of a new vector {s*(y)} whose components are 


d^{v} 



mr.i 


where 


and 


{B*(y; } 

NZxl 


(s*(y) } 
mxl 


{sj(y) } 
NZxl 


{s*(y)} 

NZxl 


{s*(y) } 


NZxl 


{ 0 } 


m-2NZ 


lie 


1 IM 


6h G 

X 


IM 

6h G 

X 


( 2 . 88 ) 


once again by following the procedure described in section 2.3, the 
equations (2.87) are reduced to a set of first order ordinary differ 
ential equations cud the new set of equations are. 
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Definitions of vectors V, S(y) and matrix A 2 are given by equations 
(2.46), (2.47) and (2.48) respectively. The vector S*(y) is a new 
addition and defined as 


S*(y) ■ 

2nxl 

The equation (2.89) can not be solved by the direct application of 
the method of recurrence relation described in section 2.3.1. It 
will require some small modification due to presence of vector S*(y). 

Let the y lines be divided into n Intervals with h^ being the 
spacing between the (l-l)st and the ith nodal points. Then, 


(s*(y) } 


fflXl 



'^1 *1-1 '’1-1 





n-i 


(2.90) 


The above equation can bo rearranged as 


-1 


[I - y A^l 


'I * T *i-i> ''1-1 


-1 




-L, 
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where 

h -1 h 

Li “ [I - T ^ T 

h -1 h 

Mi - [I - y Ail [Si - Si_i + ^ (S* + S*_i)] (2.91) 

we see an appearance of new terms in the Mi vector. Once these modi- 
fied Mi vectors are calculated, the procedure described in the section 


(2.3.1) is followed. 
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CHAPTER 3 

SOLUTION OF ELASTO-PLASTIC PROBLEM USING THE 
METHOD OF LINES 

In this Chapter a three-dimensional elasto-plastlc material be- 
havior formulation for the governing field equations in Incremental 
form Is presented. The method of lines is extended to include the 
non-linear plastic strain terms. The solution procedure Is special- 
ized for the two geometries considered herein, but can easily be ex- 
tented to other geometries. 

3.1 Governing Field Equations 

The material is assumed to be isotropic and homogeneous. The 
deformations are considered to be small and quasi-static. The Von 
Hises criterion given below is used to determine the yield condition 
at each material point. 


0 

e 




+ (0 

y 



+ <0 




2 


+ 6 (a 


2 

xy 



0 

zx 


1/2 

)1 


( 3 . 1 ) 


Subsequent yield surfaces are determined by using the rule of isotropic 
hardening. Prantl-Reuss relations are employed for constitutive 
equations relating the plastic strain rates and stresses. 

In the plastic range the strains are in general not uniquely 
determined by the stresses, but depend on the whole history of loading. 
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Therefore, the load Is applied in increments as fractions of the tota’ 
loading and the equations are written in the incremental form. Let 
be the total strain after the kth loading step. Following the conven- 
tional tensor notation it can be written as 


e 


ij 




+ 



l,j - 1.2,3 


(3.2) 


where is the strain due to the elastic deformations, t?. is the 

ij j 

accumulated plastic strain through the (k-l)th loading step and 
is the plastic strain due to plastic flow in the kth Increment. The 
equation (3.2) is rearranged as 





(3.3) 


The elastic strains are related to stresses in the following manner. 




E 


The incompressibility condition in tensor notation is 


(3. A) 


and 




0 


(3.3 


(use of the double subscript represents the summation over that sub- 
script . ) 

on making use of the equation (3.4) and (3.5), the equation (3.2) can 
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ij 






) e - ( — 5 — ) (eP + ArP 

1 - IV kk '■1 + v'' ■ ij ^^ij) 


(3.6) 


For small deformations, strain-displacement relations are. 


"iJ • s''*!.! * 


(3.7) 


Introducing equation (3.7) into equation (3.6) leads to 


^ij “ \.k 




(3.8) 


The equation (3.8) represents the constitutive relations for the elaso- 
plastic materials. The equilibrium equations for zero body force are 


0 


ij.i 


0 


(3.9) 


A set of equations similar to the Navier-Ceuc’ny equations for the elasto- 
plastic material are obtained on substituting equation (3.8) into the 
equations (3.9). These equations are. 


1 

""j.ii 1 -'2v "k.kj 



+ Ac 


P 

Ij.i 


) 


(3.10) 
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The equation (3.10) is essentially similar to the equation (2.4) ex- 
cept the addition of the plastic strain terms on the right hand side. 

For convenience the following non-dimensional quantities are Introduced, 


where 



E 

e - -22 
^xy 


£ 

E . ^ 

yz 


e 

zx 

*^zx e 

o 


(3.11) 


(3.12) 


E is the Young's modulus and is the yield stress of the material. 


°xx 

0 

0 

. 3: 

°0 

xy 

0 

o 

"hi 

0 

0 

- — 2^ 

°o 

yz 

°o 

*^zz 

^zx 

^zx 

^o 

' °0 


(3.13) 
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(3.U) 


where B Is the thickness of the crack specimen. The substitution of 
non-dimensional quantities '*1 into (3.8) and (3.10) leads 

to the following non-dimensional equations. 


(1 + V) a - 


. 3u . 9v . 3w. 

5; + 55 * 


-(r4— )Ce^ +^P ) (3.15) 

r + V ‘^XX XX 


i [vIH + (1 

(1 + V) (1 - 2 v)^ 3 x ^ ^ 


V) + v|»] 
93’ 9z 


Tt-ttt ) (3*16) 

(1 f v) ' yy yy 


°zz “ (1 + v) (1 - 2v)^ ^9x'‘' ^ 9z^ 


+ Ae^ ) (3.17) 

(1 + v) zz zz 


1 r5u . 9v, 
2(1 + v) 9y 9x^ 


— r(e^ + Ae^ ) 

(1 + v) xy xy 


(3.18) 


1 r9~ + 

2(1 + v)^!f dr 


— -i — r (e** + Ac^ ) 
(1 + v) yz yz 


(3.19) 
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zx 


1 .9w.3ui 1 ,-p .A“P\ 

i(l + Fz^ ■ (1 + V) ^^zx ^^zx^ 


and 


.2- .2- -2- 

Lii + c (— + ~) 
-2 ^ ^1' -2 -2' 
9x^ ^ 3y^ 3z 


4 5x ^3y ^ 3z^ 


+ 2C,[4(^P + ) + 4 (e** + ) 

1 3R XX XX 3y yx yx 


+ ^ (^„ + 1 


ZX zx 


+ c,(^ + ^) - c, 4l^ + §J 


_-2 .-2 
3y 3x 3z 


'2 3y‘3x 3z' 


f- 2C,[|= (fP + Ag’’ ) + « (e^ + AgP ) 

1 3x xy xy 3y yy yy 


yy yy 


■f ^ (e^ + AgP ) ] 
32 zy zy 


^ + c + ~) - c i- + ^1 
322 " %-2> S32f3x"3y-1 


+ 2C. [~ (e^ + Ac^ ) + J= (c^ + Ag^ ) 
1 3x xz xz 3y yz yz 


+ ^ (gP, + AeP )] 

3z zz zz 


(3.20) 


(3.21) 


(3.22) 


(3.23) 
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where 


C 


1 


1 - 2v 
2(1 - V) 


and 


C 


1 



V) 


3.2 Solution of Elasto-Plastic Problem 

The equations (3.21), (3.22) and (3.23) are compared with the 
Navier-Cauchy equations (2.5) through (2.7) described in Chapter 2. 

In the new equations, the field quantities and are used, 

which are now non-dimensional quantities defined by the equations 
(3.11) and (3.12). Another difference between these two sets of 
equations, is the presence of additional terms containing plastic 
strains. These differences do not alter the solution procedure 
described in the preceding chapter. The line method is still used 
to solve the equations (3.21) through (3.23). 

The equations (3.18) through (3.20) give the value of shearing 
stresses. Since the plastic shearing strains are always zero on the 
free surfaces, e.'cept in the case when the shearing traction is 
prescribed on the boundary, the plastic strain terms can be dropped in 
these equations for the boundary shearing tractions. The new equations 
will be in terms of non-dimensional quantities, but the type of 
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equations remains the same as used in Chapter 2 to eliminate the 

imaginary lines at the boundaries. This result offers considerable 

help in extending the method described in the previous chapter for 

the elasto-plastic problem. To achieve the complete equivalence we 

define the non-dimensional spacings h , h and h as 

X y z 


h 

X 



h 

y 



h 

z 



(3.24) 


where h , h and h are the line spacings for the x, y, and z sets of 
X y z 

lines as shown in Figure 3 and B is the thickness of the cracked speci- 
men. 

Ordinary differential equations for the set of x- directional 
lines are obtained by replacing x^, u^ and h^ (i = 1,2,3) by the 
corresponding non-dimensional quantities x^, u^, h^ (i “ 1,2,3) in 
equation (2.22). Hsing the matrix notations, the new equations for 
the elasto-plastic case are 




(3.25) 


^xl 


ixl fxl 


?,x2, Hxl 


ixl 
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Matrix [ic^l and vectors {u}, (r(x)} are obtained from the equations 
(2.23) and (2.2A) by using the corresponding non-dimensional quantities. 
The vectors {r^ (x) } and (r^Cx) } are new additions to the equations 
containing terms of plastic strains. They are considered functions of 
X along the x- directional lines aud expressed as follows. 


{rj(x)} ■ 
Ixl 


and 


{rP(x)} 

£xl 


2C, (r + ) 

1 XX XX 


2C, (£P + A£P ) 
1 XX XX 


2C, (eP + AeP ) 
1 XX — 


XX 


a-1 


2C. (£P + AfP ) 

1 XX XX 


2C (eP + AeP ) + ^ (eP + A£p )] 

1 9y yx yx 9z zx zx"^ ^ 

2C (eP + AeP ) + a (eP + AeP ) ] 

1 dy yx yx 9z zx zx^ ^ 




(3.26) 


(3.27) 
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The equation (3.25) is modified by following the procedure described 
in section (2.3) to obtain a set of first order ordinary differential 
equations. These equations in matrix form are 


dx 


A^U + 


d R(x) 

dJ! 


+ 


dR^(x) 

dx 



(3.28) 


The vectors R^(x) and ^^(x) are defined as 


rJ(H) 




and 


R^(5c) 


fo} 


Jcxl 




(3.29) 


(3.30) 


The definitions of other matrices and vectors used in equation (3.28) 
are the same as given by the equations (2. 42) through (2.44) but in 
this case non-dimensional quantities are employed to formulate the 
resnective matrices. 

For the set of y- directional lines. The equation (2.32) is ex- 
panded to include the plastic strain terms. Using the matrix notation 
and the non-dimensional quantities, the new second order ordinary 
differential equations for the elasto-plastic case are. 
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dM • IK. Hv) + ± (gP(Y)} + (sP(^)} 

dy2 y d9 dy 1 2 


(3.31) 


nod mXT& ourm mxl 


The matrix [K ] and vectora (v), (g(y)> are obtained by following a 

y 

similar procedure used for equations (2.33) and (2.34). The vectors 
(Pj(y)) and (B^(y)} are new additions containing the terms of plastic 
strains. They are considered functions of y along the y- directional 
lines and expressed as 


2C Al^ + AfP ) 

1 yy yy 




2C,(c*’ + Ac’’ ) 

1 yy yy 


2C- (e'’ + Ac^ ) 

1 yy yy 


v3.32) 


2CAl^ + A?^ ) 

1 yy yy 


2C,l«(eP + AeP ) -t ^ (cP^ AeP )] 

1 3x XV xy di zy ty ^ 


2C,l^(eP + A£p ) + ^ (eP + AcP )] 

1 851 xy xy 91 ty ly 


(B^(y)} - 


^ * rr <‘^xy * 


? /?P ^ A?P 


(3.33) 


2C. l^(cP + AeP ) -f ^ (fP + AeP 
1 93? xy xy df zy Sy 
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A set of first order ordinary differential equations are derived from 
the equation (3.31) by following the procedure described in section 
(2.3). Introducing the matrix notation they are written as. 


dy 


A^V + 


dy 


dS^(y) 

dy 


+ s^vy) 


(3.34) 


The vectors S^(y) 


and 


S2<y) 


are defined as 


s^(y) 


(o} 


mxl 


ISjCy)) 


mxl 


and 


sP(y) = 


(o) . 

mxl 


(3.35) 


(3.36) 


The definitions of the. othe. matrices and vectors used, in equation 
(3.34) are the same as given by equations (2.46) through (2.48), but 
in this case non-dimensional quantities are employed to formulate 
the respective matrices. 

Similarly the ordinary differential equations for the set of z- 
directional lines are obtained by extending the equation (2.38) to in- 
clude the plastic strain terms. Using the matrix notations, the new 
equations for the elasto-plastic case are 

2 “ 

- [K^l W) + i {t(5)) + i {tju)) + (3.37) 

dz 

nxl nxn nxl nxl nxl nxl 
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The matrix [K 1 is formed by using the non-dimensional spacings h 
1 # ^ 

and h in place of h and h in various submatrices of matrix [K_] 
given by the equation (2.39). The vectors (w) and {t(s)} are formed 
by using the non-dimensional quantities in the equation (2.40). New 
vectors {t^(r)> and {^^(z)} are considered functions of z along the 
set of z- directional lines. They are expressed as follows. 


ItJCi)) 

nxl 


and 


{t^(z>} 

nxl 


2C, (e^ + ) 

1 zz zz 


2C, ) 

1 zz . zz 


2C, (e^ + ) | 

1 zz zz ' 


2C-(eP + Ae^ ) 
1 zz zz 


n-1 


2C.[^(e^ + Ae^ ) + + Ae^ )] 

1 xz xz ay yz 1 

2C-[^(eP + AeP ) + |=(iP + Ae^ )] 

1 9x xz xz 9y yz yz 


2C + Ae^ + ^(e^ + Ae^ )] 

xz xz 8y yz yz 


2 C,[~(e^ + AeP ) + ^(eP + Ac^ )] 
1 3x xz xz 9y yz yz 


(3.38) 


(3.39) 


J 
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A set of the first order ordinary differential equations are derived 
from the equation (3.37) usinj’ the procedure described in the section 
(2.3). The equations in matrix form are. 


dz 


A^W + 


df(i) 

dz 


+ 


dT^(i) 

dz 


+ tP(^) 


(3.40) 


The vectors T^(z) and T 2 (z) are defined as. 


and 


Tj(i) 


T^(z) - 




(o) 


nxl 


^nxl 


(3.41) 


(3.42) 


The definitions of the other niatrioes and vectors used in equation 
(3.40) arc the same as given by equation (2.50) through (2.52), but 
in this case non-dimensional quantities are employed to formulate the 
respective matrices. 

The equations (3.28), (3.34) and (3.40) represent a set of two 
point boundary value problems, in which half of the boundary conditions 
are prescribed on one face and the remaining boundary conditions, in 
general, are given on the second face. They are either in the form 
of prescribed displacements or normal tractions. For example, for the 
set of X- directional lines, we have prescribed surface tractions on 
the face IV and face 1, as shown in Figure 4. To obtain the starting 



93 


ORIGINAL PAGE fS 
OF POOR QUALITY 

vectors for the recurrence relations method, the stress displacement 

relations given by equations (3.15) through (3.20) are used, while the 

procedure essentially remains the same as described In section 2.3.2. 

The augusented recurrence relations given by equation (2.91) are used. 

The equations (3.28), (3.34) and (3.40) are coupled differentia) 

equations. The coupling terms on the right hand side Involve u, v, w, 

their derivatives end plastic strains (l,j - 1,2,3). The solution 

Is obtained In an Iterative manner and the same Iteration Is used to 

evaluate both the dlsplacemfc.nt terms and the plastic strain terms. The 

solution of these equations yield a field solution for the displacements 

u. 's and the total strains e , e and c . The total shearing strains 
1 XX yy zz 

are obtained by using the equation 

^ij ’ 5 ^"i,j “j,i^ 


The actual evaluation will involve the use of three-point finite 
difference equations. To calculate the plastic strains we make use of 
the following equations (26). 


Ae 


Ae*’ - (2 e ’ - E ' - 

X 3c ^ X V 


S 


et 

Af 




- fx> 


AE 


Ae 


-c' - 




X 
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Ac , 

^ c 
xy Cgj. xy 


Ac 

AcP . ^ f’ 

>2 


Ac 




zx 


AS 

P r’ 

r" ^ 

et 


zx 


(3.44) 


where 


ij IJ 


n 

- E EP 


£-1 


Ij 


i.j » 1,2,3 


(3.45) 


k refers to the current increment number. 


A 


, 2 


. 2 


, 2 


e ^ [(c - e ) + (c - e ) + (e 

et 3 X y y z z 


"x> 


• 2 ,2 , 2 ^ 

+ 6(e ) + 6(e ) + 6(c ) ] 2 

xy yz zx 


(3.46) 


and 


Ac 


get - 3^^ °e,k-l 

^ 1 + 2(L±^)(_^) 

3 E Me ^k-1 
P 


(3.47) 


. do 

1 6 

The term - {- — ), , is determined from the stress-strain diagram of 
E dCp k~l 

the material. The non-dimensional effective stress o is defined as 

e 


0 - 4;t(o - 0 )^ + (0 - 0 )^ + (0 - o )^ 

e /2 X y y z z x 

-2 -2 -2 ^ 

+ 6 (o +0 +0 )1 2 

xy yz zx 


(3.48) 
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The method, applied to the present problem is known as the method of 
successive elastic solutions (26) . At first an elastic solution Is 
obtained. It Is scaled up such that the highest stressed node point 
in the grid reaches the Incipient yield condition. For the subsequent 
loading the loading path Is divided Into a number of sufficiently small 
increments. The iterative procedure for determining the incremental 
plastic strains for each load Increment is as follows, 

1) The applied load is calculated by adding the loading increment 
to the current load level. 

2) Initial values of plastic strain Increments are all assumed to 
be zero in the beginning of the load Incrertent. 

3) Obtain the field solution for the total strains ej^j(l,j ■ 1,2,3) 
by solving the equations (3.28), (3. 34) and (3.40). 

4) Calculate the modified total strains from equation (3.45) and 
evaluate the equivalent modified total strain from equation 

(3.46) . 

5) ^Ind the equivalent pl«»tlc strain increment Ae^ from equation 

(3.47) in which a . . is the dimensionless value of the 

C f i- 

equivalent stress at the end of the previous increment of load- 
ing. For the first load Increment and also for the case where 
there was no plastic flow under previous loading, o , ^ Is 
equal to the dimensionless yield stress o^, l.e., unity. 

6) Calculate r<iw set of incremental plastic strains from equations 


(3.44). 
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7) Check for the convergence of the incremental plastic strains. 
The convergence criterion used is 


-n-*-! 

yy,m 


N 


< EPS 


(3.49) 


N represents the total number of yielded nodes, while n is iteration 
n amber and m is a grid point identifier. The value of the EPS is chosen 
depending upon the desired accuracy of the solution. Repeat steps 3 
to 6 until the convergence is achieved on the incr ‘mental plastic 
strains . 

8) Sura the plastic strain increments and return to step 1. 

Once the successive approximation procedure has converged, the 
stresses are calculated at all the grid points. The schematic repre- 
seitation of the method is given below. 
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CHAPTER A 

STRESS INTENSITY FACTORS AND 
J-INTEGRAL DETERMINATION 

The stress Intensity factor (SIF) and the path independent 
J-lntegral propos^a by Rice (27) are mof.t connnonly used fracture 
mechanics parameters. Considerable effort (28,29,30) has been de- 
v'^ted to develop the techniques for an accurate determination of 
these parameters. 

In the next section a commonly used technique to evaluate the 
stress intensity factors is discussed. Finally, a procedure is 
given to determine the vao Independent J-lntegral. 

4.1 Determination of Stre.is Intensity Factors 

The stress intensity factor (SIF) is defined as 

K - Urn o (27rR)*' (4.1) 

R-K) ^ 

where 0^ is the stress ahead of ci k tip, R is the distance measured 
from the crack tip. n is the singularity of the stress field in the 
neighborhood of the crack. It was foun<i however, that due to the 
coarseness of the grid used, the usual extrapolation techniques did 
not Yield accurate results. The precise crack tip location is also 
unknown, except that it is approximately midway between two lines. 

To overcome these problems, a procedure described in reference (22) 
is followed. In this method two terms in the displacement series 
expansions around the crack tip are retained rather than one. This 
also permits to determine the actual crack tip location from the 
computed results. 
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Ar + r)3] 


(A. 2) 


where a is a function of Poisson's ratio, and r is the crack edge 

position correction, measured from the originally assumed mid point 

position. Using displacement data from the adjacent nodes to the 

Lt 

crack edge in equation (A. 2), values of the oK , ~ and r are calcu- 

i K-t 

lated for each value of z. The distance R is measured from the halfway 
point between nodes specifying boundary stresses and displacements, 
respectively. A plane strain condition is assumed to exist all through 
the thickness of the specimen, except at the surface. Correspondingly 
a would be equal to 3.56 for the plane strain case and A.O for the 
plane stress case. 

In the case of the curved crack fronts, the v- displacement values 
are extrapolated along the normal directions tr the crack front, with 
the help of nodal displacements. These values are used in equation 
(A. 2) to evaluate the local stress intensity factors through the thick- 
ness of a cracked specimen. 

A. 2 J-integral Determination 

Rice (27) developed a path independent integral, known as the 
.-integral in the literature. This integral is associated with the 
change in energy of a body due to the crack growth and it is expressed 
as 


J 


3u 

(W(e) dy - T^ ^ ds) 


(A. 3) 
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Here r Is a curve surrounding the notch tip as shown in Figure 9(a). 
The contour integral is evaluated in a counter clockwise sense, 
starting from the crack surface. The strain energy density W(c) is 
defined as 

£ 

mr. 


mn 




(4.4) 


o 


and the traction vector is 


T 


i 


"j 


(4.5) 


u^ is the displacement vector. 

The consideration of singular terms associated with the stress environ- 
ment near the crack tip in a linearly elastic body, in the evaluation 
of equation (4.3) leads to well known relationships as obtained by 
Rice (27). These relations are. 


J 

J 


1 - „2 
E 


E 


(for plane strain) 
(for plane stress) 


(4.6) 


The equations (4.6) could be used to determine stress Intensity fac- 
tors by computing J-values which can be obtained without a detailed 
knowledge of the stress and strain field, very near the crack tip. 

I By taking advantage of geometric and loading symmetry about the x 

^ axis, retangular paths are chosen to calculate J- integral values. 

[ For a rectangular path shown in Figure 9(b) the contour integral for 
J can be written as, 
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i: 


3v 


J - 2 1 [W(e) - “ o,y 3 ^ Idy 


+ 2 1 (a e„ + 1 ^) dx 


+ 2 


xy X y 9x 


3v, 


W(e) - o t - a ^]dy 
XX 3x 


where 


W(e) 



de +0 
X y 


de + a 

y 2 


de 

z 


+ 2 o de + 2 o de + 2 o 
xy xy yr yz zx 


(4.7) 


de ) (4.8) 

zx 


The coordinates of the points shown In Figure 9(b) are used as Inte- 
gration limits in equation (4.7) and the integrations are calculated 
using the trapezoidal formula of numerical Integration. 
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A computer program was written to solve the elastic problem for 
the different cracked specimens. This program was later modified to 
Include the plasticity terms. Since the ordinary differential equa- 
tions are coupled, a successive approximation procedure was used to 
obtain the solution. The computation for ail the examples were per- 
formed on IBM 370/3033 time sharing computer. 

5.1 Single Edge-Notched Tensile Specimen 

The new formulation of the MOL was first applied to a single edge 
notched specimen in tensile loading. The specimen shown in Figure 2 
has dimensions, W/a = 2, B/a = 3, L/a = 3.5, where a, B, W and L are 
the specimen crack length, thickness, width and length rcspectivcl v . 

Figure 10 shows a plot of the <iimenslonless maximum crack opening 
versus the number of grid lines. We note the following: 

1) This displacement approaches smoothly an asymptotic value. 2) The 
value obtained in reference (31) using the old formulation of MOL, is 
in error by approximately 17% whereas the errcir using the same number 
of lines with the new formulation is onlv about 3 percent from the 
asymptotic value. Furthermore, 12 x- g.rid lines were the most that 
could be used in reference (31), before the solution diverged. 3) A 
comparison with the results obtained in reference (32) using finitti 
elements, indicates that the present solution using 16-x lines, is 
essentially correct since it was indicated in reference (32) that the 
value obtained there may be as much as 25 percent low. 
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The dimensionless stress Intensity factor for this specimen at 
the mid thickness is shown in Figure 11 versus the number of x grid 
lines. Again, the value obtained in reference (31) using the old 
formulation is in error by about 8 percent, whereas the new formulation 
with the same number of grid lines is in error by about 3 percent from 
the asymptotic value. 

Figure 12 shows the stress distribution ahead of the crack tip 
at the mid thickness for different grids. Note that the same distri- 
bution is obtained with 12, 14 or 16 x grid lines. The larger the 
number of x grid lines, however, the closer the crack tip is approached, 
and the more accurate is the stress intensity factor calculation. 

The variation of the stress intensity factor through the thickness 
of the specimen is shown in Figure 13. The SIF decreases approximately 
15 percent in going from the center of the specimen to the surface. 

5.2 Compact Tension Specimen 

The standard compact tension specimen (CTS) shown in Figure 1(a) 
is considered next. The dimensions of this specimen are W/a » 2, 
b/a * 1, L/a = 2.4. The load was assumed to be parabolic shear load 
."long the pin load line as shown in Figure 1(b). 

Figure 14 shows the dimensionless crack opening displacements for 
this specimen. As is seen, there is very little difference between the 
midplane and surface displacements. The experimental value from refe- 
rence (25) is also shown and it is seen that good agreement is obtained. 

The variation of stress Intensity factor through the thickness of 
the compact tension specimen is shown in Figure 15. Also, shovm in 
this figure are some of the results obtained by finite elements in 
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reference (11), which treats the Battelle Benchmark geooetrls, (33). 
Several different methods were investigated In reference (11) for 
e**aluating the stress Intensity factors from the finite element results. 
These were d^.splacement substitution, modified displacement substi- 
tution, virtual crack extension and contour Integration methods, desig- 
nated as methods D, E, V, and J respectively in reference (!?' 

In the displacement substitution method, the displacement at the 
vertex and mid-edge nodes In the special elements on the face of the 
crack were used. Since a coarse mesh was used In reference (11), an 
extrapolation was done on the values of SIF obtained by displacement 
substitution from the vertex and mid-edge nodes and the method Is 
named as modified displacement method. 

In the virtual crack extension method the evaluation c f the energy 
change is made corresponding to a small adjustment in the position of 
any point on the tip of the crack In any direction. As shown by Rice 
(27), the component equals the rate of change of energy per unit 
area of crack extension at a point in the direction of the x- axis. 

The path Independent J- integral proposed by Rice (27) Is used in the 
method J. 

The results of methods D, E and V have been plotted In Figure 15 
together with the results obtained herein using the MOL. The results 
obtained from MOL lies approximately midway between the methods V and 
E. The result from the two-dimensional plane strain solution given In 
ASTM E399-78(5) Is also shown in Figure 15. The MOL result is only 
0.6 percent higher than this result at the canter where plane strain Is 
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expected to exist. The two dimensional solution Is obtained by boun- 
dary collocation method and believed to have an accuracy of 0.2 percent 
fox a CT specimen. Note also the drop of approximately 15 percent In 
SIF value as one moves from center to the surface of the specimen. 

From the above results It Is apparent that the line method of 
analysis as modified herein appreciably Increases the accuracy and 
reduces convergence problems . 

5.3 Curved Crack Front Specimens 

The method of lines was next applied to curved crack front (CCF) 
problems. The compact tension specimen In tensile loading Is consi- 
dered first. An accurate geometrical description of a curved crack 
front Is net possible If a coarse grid such such as NX * 12 » NY ■ 6, 

NZ ■ 5 Is used. Therefore, two fine grids of NX “ 20, NY ■ 7, NZ • 5 
lines and NX* 20, NY *7, NZ*7 lines were used In the subsequent 
analyses. Since approximate parabolic curves were used In the study It 
was difficult to assign a single measurement for the crack front curva- 
ture. To describe this curvature, therefore a parameter called crack 
tunnel depth is defined as the difference between the crack lengths at 
the center and the surface of the specimen. Thus an Increasing crack 
tuunel depth reflects an increasing crack curvature. 

Figure 16 shows a plot of the local stress intensity factor for 
a CT specimen In tensile loading. Curve 1 shows the SIF variation 
with thickness for a straight crack front. For this case maximum SIF 
occurs at the center of the specimen and it drops by 13 percent on 
the surface. This trend Is in accordance with the previously reported 
results of Raju (10). For the crack tunnel depth of 0.1 it Is seen 
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that SIF decreases at the center by 9.5 percent, while at the surface 
its value increases by 20.5 percent. Similar type of variation is 
reported by Pereira et. al. (14). However, a direct comparison of 
the results is not possible because they used a different crack and 
specimen size. For higher crack tunnel depths the maximum SIF value 
no longer occurs at the surface. It shifts to the adjacent interior 
point of the grid. This trend is not reported in reference (14), 
probably because a crack length (a/W) of 0.25 was used which is 1/2 
of the one used in the present study and more importantly three layers 
of elements were used in the thickness direction, making it probable 
that the maximum SIF value was missed. 

Figure 17 shows the variation of non-dimensional SIF for the CT 
specimen with assumed parabolic loading along the crack line instead 
of tensile loading. The dimensions of the specimen are the same as 
before. Once again for the crack tunnel depth of 0.097 the local SIF 
value decreases at the center while its value rises at the surface of 
the specimen. For the higher crack tunnel depths the trend of lower- 
ing the SIF value at the center of the specimen continues. The 
maximum SIF no longer occurs at the surface, but it shifts in the 
interior of the specimen. This is similar to the observed trend for 
the tensile loading. 

A possible explanation for these trends could be given by consi- 
dering the various effects wliich have direct bearing on the SIF values. 
The rise in the SIF value at the surface is mainly due to the presence 
of a crack front with curvature. However, there are two other effects 
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which are responsible for the reduction in the SIF value at the surface. 
Firstly, the existence of plane stress condition at the surface leads 
to a drop in SIF value (an approximate drop of 15 percent was observed 
for a straight crack front). Secondly, due to the presence of a curved 
crick front, the crack length decreases at the surface as compared to 
th3 center, leading to a reduction in the SIF value. As it could be 
sein that for the smaller crack tunnel depth the effect of curvature 
is dominant over other effects and thus the maximum SIF value occurs at 
the surface. However, for the curved crack fronts with larger crack 
tunnel depths, the effect of reduction in the crack length seems to be 
dominating. Consequently the SIF value at the surface no longer remains 
maximum. A combination of the three effects, mentioned earlier, deter- 
mi les the location of maximum SIF value through the thickness of the CT 
sp icimen. 

Figure 18 shows the variation of center, surface and thickness 
av irage stress intensity factors with increasing crack tunnel depth for 
a "T specimen. The stress Intensity factors are non-dimensionalized 
with the respective values of SIF occuring for a straight crack. The 
thickness average stress intensity factor is calculated by taking 

average of the SIF values across the thickness. It is reflective of 
the average stress environment existing near the crack tip. Note that 
non-dimensional SIF value at the center of the specimen keeps on de- 
creasing as the crack tunnel depth increases. However, the value of 
SIF at the surface Increases only up to a certain crack tunnel depth. 
The thickness average SIF also decreases with increasing crack tunnel 
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di'Pth. Similar behavior for the tliickness average SIF was observed 
by McGowan (12) in his study on a single edge notch specimen. This 
clearly shows that the presence of curved crack front can si<^nif icantly 
change the alue of SIF. The lowering In the value of SIF indicates a 
reduction in stress concentration near the crack tip. Consequently, the 
specimen can withstand higher loads before showing any appreciable crack growth. 

In Figure 19 the change in the thickness average SIF is plotted 
against a non-dimensional curvature parameter n. The value of n corres- 
ponds to a difference between the average crack length and the 

surface trace of the crack divided by average crack length (a ). 

The calculation of a^^ is done by computing an average of the three 
crack measurements taken at the center of the crack front and midway 
between the center and the end of the crack front on each side. This 
definition of a is in accordance with the ASTM standard E399-78 (5). 

The standard restricts the length of either surface tra:e of the crack 
to more than 90 percent of the average crack length (a^^) • Based on 
this criterion, the ASTM region is marked on the Figure 19. Any test 
in which the specimen shows a surface trace within the marked region, 
will be considered a valid test. It could be seen that for the case 
when there is a surface trace of 90 percent of a^^ in the specimen, 
there is a 7.0 percent lowering in the thickness average stress intensity 
factor. This shows that the use of plane strain formula given in 
references (5) will result in an overestima'-ion of the fracture toughness 
of a material by 7.0 percent. 

On the basis of experimental results and three-dimensional analysis 
of Pereira et al. and Neale (14,15), recently, an amendment has been 
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incorporated (40) to extend the limit on the surface trace of the crack 
to 85 percent of the average crack length. This corresponds to a n 
value of 15 percent. A drop of 12.5 percent in tua value of thickness 
average SIF can be observed from the Figure 19 for this value of r\. 

This is an Increase o.. 5.5 percent over the previous value and implies 
that the use of two-dimensional plane strain formula will result in an 
overestimation of the fracture toughness by 12.5 percent. These figures 
provide a fairly good estimate of overestlmatlon of the fracture tough- 
ness if the CCF is present in the test specimen. 

5.4 Effect of Plastic Flow 

The increase in SIF values at the surface of the CT specimens with 
different curved crack fronts Indicates that the crack Initiation will 
start from the surface of the specimen. This is contradictory to the 
experimentally observed behavior of crack tunnelling. It has been 
suggested that certain amount of plastic flow takes place at the surface 
which results in lowering of the SIF values there. To verify this an 
elasto-plastlc stress analysis was carried out for the CT specimen. 

The material considered in this study was M-5083. A stress-strain 
diagram for this material is shown in Figure 20. The yield stress (o^/ 
is taken as 12.417x10 psi and the value of elasticity modulus (E) is 
estimated from the test data provided in reference (34), The plasticity 
solution was obtained by following an incremental solution procedure. 

The total numer of increments needed for complete elasto-plastlc analysis 
were between 20 and 26. The Incremental load ~ was taken approximately 
equal to 0.075. For one increment, the complete elasto-plastlc solution 
took approximately 6 minutes of CPU time on the IBM 370-3033 computer 




109 


ORIGINAL PAGE IS 
OF POOR QUALITY 


and for the convergence criterion described in Chapter 3, a fixed value 
of 0.001 was taken for EPS. 

In Figure 21 the variation of non-dimensional crack opening dis- 
placement along the width of the CT specimen ii shown. These displace- 
ments are extrapolated to determine the crack mouth opening displacements 
of the specimen. The displacements considered herein are at the center 
of the specimen for a straight crack front. The different curves re- 
present increasing applied load levels. 

An experiment was performed to obtain the load versus crack mouth 
opening displacement (CMOD) for the compact tension specimen at NASA 
Lewis Research Center. The specimen was made according to the ASTM 
Standard E399-78 (5) . It was loaded in an MTS automated testing system 
equipped with a load cell capacity of 2C Kips. The loading was stroke 
controlled having a rate of 0.02 inch/minute. The main feature of the 
testing system included a closed-loop, servo controlled hydraulic oper- 
ation. It was fully computer operated. A clip guage satisfying the 
requirements of the ASTM £399 standard was mounted on the specimen to 
measure the CMOD. 

The stroke controlled software developed at Lewis Research Center 
was used for the data acquisition. Tlie load and the corresponding CMf)D 
constitute a data point. These data points were stored on hard disc 
with an acquisition rate of 1000 data points in half a second. A real 
time plot of load versus CMOD was cirawn on a flat bed plotter. 

The load versus CMOD results of the experiment were also printed 
in a tabular form on a line printer. This numerical data is utilized 
to draw the experimental curve on Figure 22. Also shown on this figure- 
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Is the results obtained from the MOL. At Incipient load level where 
plastic flow begins, the difference between the experlnental displaceTiient 
and the calculated value la 5 percent. There seems to be a fairly good 
correlation between the results obtained by the experiment and the cal- 
culated value. 

A complete elasto-plastlc analysis was carried out for three 
different crack fronts. These ciack fronts are shown in Figure 23. 

Tie first crack front (CCFO) Is a straight crack while tlM other two 
(CCFl, CCF2) have some curvature. The growth of the plastic cone with 
load for a straight crack front is shown in Figure 2v . The growth Is 
shown at three different locations through the thickness of the rpeclmen, 
at the center, at the surface and some plane adjacent to the surface. 
Location adjacent to the surface was chosen because of maximum plastic 



conventional two-dimensional plane model. Note also the growth of the 
plastic zone at the back surface of the specln^n. In this zone the 
plastic strains are compressive due to the bending component c£ the 
load. For the two other crack fronts CCFl and CCF2 the development of 
the plastic zone growth Is shown In Figures 23 and 26. 

To compare the plastic zones for different crack fronts, the highest 
load level of 2.3 kips was chosen. At this load level the plastic zones 
at the three different locations In the thickness of the specimen are 
plotted In Figure 27, At the center, the plastic zone sizes are approxi- 
mately the same for CCFO and CCFl, while It Is considerably smaller for 
CCF2. On the surface the plastic zone Is bigger for the CCF2. Note 
also the presence of a small elastic region within the plastic zone. 



Ill 


OF POOR 


iu l-J 

h'JTY 


This growth of the plastic zone Is anticipated for CCF2, due to 

the rise in stresses at the surface ol the specimen. The plastic zon»' 
growth is smaller for the CCFO at the surface as compared to CCFl , 

The stresses are relieved due to the large plastic flow at the 
surface, when the curved crack fronts are present. This explaines the 
beginning of crack advancement at the center of the specimen. The 
presence of triaxial constraint due to the stress at the center also 
contributes to t'le tunnelling behavior of the crack. 

5.5 Calculation of J-lntegral 

The J -integral proposed by Rice (27) plays an important role in 
the non-linear fracture mechanics. It was shown that for linear and 
non-linear elastic materials, this integral is path Independent. 
Furthermore, based on energy considerations, it was proven to be equi- 
valent to energy release rate per unit crack extension. This conclusion 
is very important because it removes the need of accurately determining 
the stresses near the crack tip. J is also used in predicting the 
el. St ic-j lastic crack growth (55). In the present study J values are 
calculated for different curvi'C crack fronts to analyze L.. • behavior 
of crack growth in the comp.ict tension spec men. 

Nine differ«Mit paths were chosen to evaluate' the J values at 
different thlcknc*sscs of the CT spe'" 1 '. 11011 . Those ptiths are shown In 
Figure 28, The non-dimensional J values at three dllfcreni locations 
are plotted in Figures 2d, 30 and 31. Tlicse ’.esiilts are for .1 stra^K.ht 
crack front at increasing load levels. The following observations 
are made, 1) For the loads which are close to the load for incipient 
plastic flow, the J values remain nearly the same for different paths. 
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2) For the higher loed levels soac perturbation is observed in the 
J-values along the different paths. Sinllar trends are reported in 
reference 36 in which a t%io>dliDensional elastic-plastic analysis of a 
CT spe'iaen is carried out. 3) As we nove from center to the surface 
of the specimen the path independence property of the J-lntegral de- 
teriorates. 

The use of two-dimensional definition of the J-integral In a 
three-dimensional analysis is debatable but as reported in reference 
(37) that due to the symmetry of the specimen it is enough to consider 
only one component of the general three-dimensional vector J. An 
equation to evaluate was given, which is similar to the two-dimen- 
sional definition of the J-integral with the addition of a surface in- 
tegral term. The detaU.: of the equation can be found in Appendix D. 

It was further shown in the same study, that the surface integral does 
not significantly change the numerical value of J. Only at the surface 
does its evaluation have some effect, on the J-lntegral values. Based 
on this, in the present study the two dimensional definition of J is 
adopted. 

Figure 32 shows the different paths used to evaluate J-integral 
values for the curved crack front CCFl. All the paths are the same as 
used before with the exception of the paths chosen at the surface. 

Since the crack length is smaller at the surface due to the curvature, 
all the paths are shifted toward the left. Once again the J values are 
plotted for different locations through the thickness of the speclsmn 
in Figures 33-33. It is observed that for lower load levels the J-lntegral 
is essentially path Independent but for the higher loads this property 
does not hold good. 
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In Figure 36 different paths, used for the curved crack front 
2 are sho«m. Note the shifting of paths towards the left at locations 
g ■ 0.375 and c > 0.3 due to the presence of the curved crack front. 

The J-integral values are plotted for different paths and locations 
In Figure 37-39 for CCF2. Once again, it Is observed that path 
independent property of d-lntegral Is closely followed only for the 
loads which are closer to the lo.-id .at which plastic flow begins and 
at the center of the specimen. 

From the above results it Is cleai that J-integral has significant 
path dependence ioaediately adjacent to a crack tip under small-scale 
yielding conditions in an elastic-plastic material. Parks (38) and 
HcNeeklng (39) have also reported similar results. This may be due to 
the unloading which occurs at the nodes near the crack tip because of 
the plastic flow. These observations very seriously challenge the role 
of J as .1 p.trameter characterising the crack-tip stress field , at 
least for materials OKidelled by the Von Mises flow theory. 

In view of these observations it was difficult to choose a repre- 
sentative .I-lntegral value for each of the loading cases. A reasonable 
choice for J can be obtained by considering Figure AO. The non-dlmen- 
sli'nal values of the stress intensity factor were calculated uslnp the 
J values oht.alned from the path 1. .Similar calculations are mad*' !<o 
the SIF using the average J values which are obtained by averaging all 
the J values for path.s 3 through 9. These numerical values of SIK and 
the SIF values obtained from crack opening displaceim*nt s are plotted In 
Figure 40. It can be seen that the SIF values obtained t rora the J-lnfcgr.il 
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values for the path 1 are very close to the SIF values froai crack 
opening dlsplacetoent . While the average J values do not yield an 
accurate result as one proceeds towards the surface of the apeclaen* 
Because of this reasor) the results from the path 1 for the J-lntegral 
are accepted as good estimates for further analysis. 

The non'-dlaenslonal J-lntegral values as a function of load, ob- 
tained from the use of path 1 are shown In Figures 41-43. Curves for 
the three different crack fronts are shown In the figures. From figure 
41, It could be seen that at the center of the specimen, the J values 
for CCFl and CCF2 are lover than the values for CCFO, but at a load of 
2.3 kips the J values of CCFO and CCFl coincide while for CCF2, it 
continue to remain lower. On the surface, the J values remain higher 
for CCFl and CCF2 for the smaller load levels as shown In Figure 43. 

But, as load increases this no longer remains true and the J values 
go down considerably for CCFl. The J values for CCF2 also fall below 
the J values of CCFO. This shows that due to the plastic flow, J- 
integral value decreases at the surface for the curved crack fronts 
as compared to straight crack front. 

The variation of J obtained from path 1 Is plotted through the 
thickness on Figures 44-47. For the straight crack front the J values 
for the lover loads are maximum at the center and minimum at the surface 
as shown In Figure 44. With the Increase In load this trend Is slightly 
modified. The J-lntegral value at the center still remains the maximum 
but the minimum value shifts from the surface to some Interior point. 

For CCFl and CCF2, the J value Is maximum at the surface for the loads 
close to load at which plastic flow begins. For higher loads this trend 
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reverses for CCFl. The aaxij»m value of J shifts to the center of the 
specimen while for CCF2 the maximum values do not occur at the center 
but shifts to an interior point. 

The existence of a critical J* value is assumed at which crack 
advanceTOnt starts. From Figure 45, it is obvious that for CCFl after 
some plastic flow, the crack advancement will start at the center of 
the specimen, not at the surface as predicted by a purely elastic analysis. 
Similar behavior was deduced from the approximate elastic-plastic model, 
developed by Neale (16) for a CT specimen. This is in full accordance 
with the crack tunnelling behavior observed during the experiment. For 
CCF2 the maximum J value will occur at the center of the specimen, as 
the shift in maximum J value towards the center could be seen in Figure 
46 with the increase in load. 
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CILLER 6 

SUMMARY AND CONCLUSIONS 

An Improved formulation of the method of lines (MOL) is presented. 
The five point finite difference formulas are Introduced to achieve 
more accurate results. The resulting ordinary differential equations 
are solved by a recurrence relation method. It is a veil suited method 
for solving tvo point boundary value problems. Two specific geometries 
namely the edge notch specimen and the compact tension specimen are 
considered. For both the geometries, complete field solution for the 
stresses and strains were obtained. 

To establish the convergence characteristics of the newly Improved 
MOL, the field solutions were obtained for edge notch specimen for 
different grid sizes. On comparing the maximum crack opening and the 
stress Intensity factors for different grid sizes, it was found that 
the solution converges to an asymptotic value. Even the coarse grid 
such as NX » 12, NY > 6, NZ * 5, yielded fairly good results and the 
CPU time was only of the order of 2.5 minutes. 

For the compact tension specimen maximum crack opening displacement 
was compared with experimental results and were found to be In good 
agreement. The SIF value at the center of the compact tension specimen 
was only 0.6Z higher from the value given In the ASTM standard (reference 
(5)). 

Only smaller grid sizes such as iDC « 12, NY ■ 8, NZ ■ 8 could be 
used In the old method of lines. For the bigger grid sizes the solution 
became Instable. This places a serious restriction on the shape of the 
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crack fronts which could be used to study the effect of crack front 
curvature on the local stress intensity factors. By introducing the 
new modifications in the MOL« it became possible to use bigger grid 
sizes such as NX - 20, NY - 7, NZ - 7. Presently, only the core size 
and the CPU time seem to limit the size of the grid. 

The complete solutions were obtained for the compact tension 
specimen in tensile and shear loading, containing different crack 
fronts. For the purely elastic case, it was found that as the crack 
front curvature increases the SIF value at the center of the specimen 
decreases %*ile increasing at the surface. For higher values of crack 
front curvatures the maximum value of the SIF occur at an interior 
point located adjacent to the surface. These results indicate that 
for the specimen containing a curved crack front, the crack growth will 
Initiate at the surface of the specimen. This conclusion is in direct 
contradiction with the experimental observations in which a tunnelling 
behavior is observed, but can be explained by the presence of plastic 
flow. 

A thickness average SIF value was computed. It is assumed that 
it reflects the average stress environment near the crack edge. On 
the basis of this, it was estimated that use of the ASTM formula 
(reference (5)) will lead to an overestimation of the fracture toughness 
by 7 percent if the curved crack front present in the compact tension 
specimen Just satisfies the ASTM limit on the surface trace of the crack, 
providing no plastic flow occurs. It was further estimated that the 
proposed amendment in the ASTM standard on the surface trace of a 
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crack will lead to a maxitnun overeat Imat Ion of the fracture toughness 
by 12.5 percent. 

To Investigate the effect of plastic flow In the compact tension 
specimen, the equations of the method of lines are augmented to in- 
clude the plasticity terms. Complete elasto-plastlc analysis were 
carried out for three different cases of crack fronts, which Include 
the case of straight crack front. 

To check the accuracy of the elasto-plastic analysis, the load 
versus crack mouth opening displacements were compared with experimen- 
tal results. The experiment was conducted at NASA, Lewis Research 
Center. The two results were found in good agreement. 

The growth of the plastic zones are compared for a straight crack 
and two curved crack fronts. The maximum plastic zone occurs at the 
surface for the curved crack front with the maximum curvature. This 
qualitatively explains the reason for the initiation of crack growth 
at the center of the specimen. Due to the large plastic flow at the 
surface, the stresses are relieved. 

To further investigate this, the J-integral values are computed 
at different locations for each of the crack front. In general, it is 
observed that the J values are path independent for the lower loads. 

As the extent of plasticity increases, this property of path independence 
breaks down. It was also found that as we move from the interior 
to the surface, this property deteriorates. In view of these results, 
one specific path was chosen to compare different J values. The choice 
of the path was based on a comparison between the SIF values obta'ned 
from the J-integral values and the SIF values computed using the crack 
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opening displacements. For this path it was observed that at the load 
levels close to the elastic conditions the J values are higher at the 
surface and lower at the center for the curved crack fronts. As the 
load Increases, the trend reverses itself and for the curved crack 
front with small curvature, the maximum of J value occurs at the center. 
This provides a quantitative explanation of crack initiation at the 
center of the specimen. This conclusion can not be drawn from a purely 
elastic analysis. 

6.1 Concluding Remarks 

The new improvements in the method of lines, have considerably 
enhanced the accuracy and the stability of the method. Converged 
results can be obtained by using ;alatively coarse grids. It was 
observed that increase in the number of lines in one direction only, 
can lead to an instable solution. In some cases this instability 
could be removed by adjusting the spacings, used for the application 
of the recurrence relation method. These spacings were selected in 
such a way, that they were approximately equal to each other in all 
the three directions. This was established purely on the basis of 
numerical experimentation. At this stage, a rigorous mathematical 
analysis of error estimates and stability of the method is desirable. 

In the present work, one particular constant grid size, was cnosen 
for each direction. This use of the constant grid size to approach 
crack tip as close as possible, leads to a l?’‘ge number of equations 
to be solved. In any future work this problem can be alleviated by 
using a varying grid size. 
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Th« prM«nt formilatlon of the elastlc-pleetic problea ie besed 
on the euccessive eleetlc solution acthod. It Is observed that rate 
of the convergence slows down as the extent of plasticity Increases. 

The solution procedure start to diverge when the plastic som starts 
growing in the coapresslon sons of the compact tension specimen. This 
shows that present formulation cannot predict large non-linear effects. 
It will need some modifications to obtain solution for large plasticity 
effects. 

The Iterative scheme of the method of lines was based on the 
successive approximation procedure. This scheme la simpller to 
adopt* but needs more Iterations. To cut down on the Iterations* for 
the elastic case a successive over-relaxation (SOR) parameter was tried. 
An Improved guessed solution was obtained by cond>lnlng the solution 
for the current and previous Iteration with the help of a SOR parameter. 
This Improved guessed solution was used to carry out the next iteration. 
In the case of a grid* MX* 20* NY ■7*NZ*S, the numbpr of Iterations 
were reduced from forty-seven to twenty-eight* when a SOR parameter 
of value tvo was used. This leads to a considerable saving In the 
computer time. Such parameters could be studied In more detail to 
accelerate the convergence rate of the present method of line. 
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Figure 1(a). Standard Compact Tension Specimen. 



Figure 1(b). Idealized model of Standard Compact Tension 
Specimen. 


Figure 1. Standard Compact Tension Specimen and its 
idealized model. 
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Figure 2, Single Edge-Notched specimen with tensile loading. 
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a) An arbitrary contour in an non-linrar elastic body 
containing a crack. 



b) A rectangular path with coordinate numbering,. 
Figure 9. Path for the determination of J-integral. 
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Figure 10. Variation of maximum crack opening displacement 
with number of x-grid lines for single edge- 
notched tensile 8peclmen,NY«NZ"8. 
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Figure 11. Dimensionless SIF at the center of single 
edge-notched specimen versus number of 
x-grid lines, NY*NZ“8. 
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Figure 13. Variation of nondimensional stress intensity 
factor through the thickness for single 
edge-notched tensile specimen. 
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Dimensionless stress intensity factor, K^bJw/P 
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2D . In. ASIM E399-78 (5) 
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O Method D Ref. (11) 
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Figure 15. Variation of dimensionless SIF through the 
thickness of the compact tension specimen. 


Diniensionler?s stress intensity factor, K bJw/P 
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Figure 16(a). Different crack length measurements. 



Distance along the thickness, z/B 


Figure 16. Variation of d imens ion i ess SIF wltii thickness 
for different crack tumel deptlis for compact 
tension specimen under tensile loading. 
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Figure 18. Variation of center, surface and thickness average 
stress intensity factor with increasing crack 
tunnel depth. 
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Figure 19. Percent variation in average stress intensity 

factors between straight crack and curved crack 
fronts . 
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gure 21. Variation of non-dimensional crack opening displacement 
along the width of the compact tension specimen. 
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Figure 22. Load versus maximum crack opening displacement 
for the compact tension specimen. 
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Figure 25. Growth of plastic zone with Increasing load 
for curved crack front 1 at three locations 
of the CT speclnen. 
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Figure 26. Growth of plastic zo-ie with Increasing load 
for curved crack front 2 at three locations 
of the CT specimen. 
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Figure 27. Plastic zone enclaves for different curved 
crack fronts at three locations of the CT 
speciaen (P«2.5 Kips). 
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Figure 29. Non-dlaensiona' J-lntcgral values for dlffarant 
paths at the center of the CT speclnen with a 
straight crack front for increasing load. 
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Figure 30. Non-dimensional J-integral values for different 
paths at midway between center and surface of 
the CT specimen with a straight crack front for 
increasing load. 
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Figure 31. Mon<>diacn8icmal J-integral values for different 
paths at the surface of the CT specimen with 
straight crack front for increasing load. 
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Figure 32(a). Different paths for the plane located 
at 2“0.0 through 0,375 



Figure 32(b), Different paths for the plane located 
at z*0,5 


Figure 32. Different paths used for evaluating J-integral 
for curved crack front 1 




Dlnensionless J-lntegral, EJ/o^B 




Figure 34. Non-dimensional J-integral values for different 
paths at midway between center and surface of 
the CT specimen with curved crack front 1 for 
increasing load. 
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Figure 35. Non-diaenelonel J-integral values for different 
paths St the surface of the CT speciinen with 
curved crack front 1 for increasing load. 
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Figure 36(a). Different paths for the planes located 
at z«0.0 through 0.25 



Figure 36(b). Different paths for the plane located 
at 2=0.375 



Figure 36(C). Different paths for the plane located 
at 2=0.5 


Figure 36. Different paths used for evaluating J-integral 
for curved crack front 2 
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Figure 39. Non-dimensional J-integral values for diff( 
paths at the surface of the CT specir.'in wi 
curved crack front 2 for increasing load. 
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Figure 40. Variation of dimensionless stress intensity 
factor through the thickness of the compact 
tension specimen . 
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Figure 41. Non-dimensional values of J-lntegral at the 
center of the CT specimen for Path 1 versus 
the applied load for three different crack 
fronts, s-0.0 
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Figure A2. Non-dlnensional values of J-integral at aidvay 
between center and surface of the CT specioen 
for Path 1 versus the applied load for three 
different crack fronts, z-0.25 
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Figure 43. Non-dimensional values of J-integral at the 
surface of a CT • . - cimcn for Path 1 versus 
the applied load for three different crack 
fronts, *"0.5 
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Figure 44. Variation of non-dimensional J-integral values 
through the thickness of a CT specimen with 
a straight crack with increasing appled load . 
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Figure 46. Variation of non-dimensional J-integral v .• s 
through the thickness of the CT specim h 
curved crack front 2 with increasing applied load* 
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APPENDIX A 


Formulation of the Governing Ordinary Differential Equations 
for the Lines Located at the Surface and 
Adjacent to the Surface of the Cracked Specimen 

For boundary surface lines and lines adjacent to the boun- 
dary surface lines, the difference expressions for the second 
derivative will involve imaginary lines outside the boundary. 
Since three-dimensional elasticity problems have three boundary 
conditions at every point of the bounding surface and a second 
order ordinary differential equation needs only two conditions, 
the shear stress at the boundaries is used to eliminate the imag- 
inary lines outside the surface while the condition of the pre- 
scribed normal traction or displacement will be enforced through 
the constants of the homogeneous solutions. 

A.l Derivation of Ordinary Differential Equation for x - Direc- 
tional Lines 

Let us take the x-dlrectlonal line which is formed by the 
Interaction of x - r and x - y coordinate planes of a solid (Fig- 
ure 6). The ear stress for ♦’he free surface is given by 

0 

yx 


x -2 Coordinate plane ■ 0 


(A.l) 


or 3u/9y + 3v/9x » 0 


(A.2) 
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At line 1 using the iToaglnary displacement of the ficti- 
tious line 1^'' the equation (A. 2) can be written as 

^2 ~ ^1 ^ 1 . dv I 

2 h 5>x ' 1 dx '1 

y 

or Uj^^^ • U 2 + 2 hy ( dv/dx) | (A. 3) 

The plane x-y Is one of the symmetry planes and u displace- 
ments are isymmetrlc with respect to this plane. 

f V 

The fictitious displacement ^ Is used to write a forward 
finite difference equation for the partial derivative with respect 
to y, while two fictitious displacements ^ and u^ ^ are used to 
v;rlte a central difference equation for the partial derivative with 
'■espect to z. The ordinary differential equation for line 1 is 
written as 

( d^u, /dx^ ) + ( C, /12 ) [11 u/^ - 20 u, +6 u. + 

11 y 1 1 L 

4 U 3 - u^ ] + ( C^/12 ) [- + 16 

"30u, +16ujjY + i- U 2 NY*ti] 

S. ~hr ''^i 


Due to the s>nmnetry considerations, 


. flz 

“1 " “2 NY + 1 


. f 2 z 

^ “ Vtl 


(A. 5) 
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on tubstituting nquatlont (A. 3) and (A. 5) Into aquation (A. 4), wa 
obtain 

( d^Uj^/dx^) +( €2^/12 h^) [- 20 Uj^ + 17 U2 + 4 Uj - u^] 

+(€2^/12 + + 

- C2 -3— (V + w)^ -(11 Cj/6 hy) (dv/dr)| ^ CA,6) 

Slallarly for lina 2, by making use of tha flctltloxia line 
fy 

1 ^ and aynnwtry tha ordinary differential aquation can be written aa 

(d^U2/dJf^) + (€2^/12 hy) [ - + 16 U2^ - 30 U2 + 

16 U3 - u^3 +(C^/12 h^> C-30 u2 + 32u2^^2" 

2 “2 NY + 23 • ^2 4;^ ''>2 

Note tha use of tha central difference equation for the par- 
tial derivative with raapact to y. The fictitious displacement 
fy 

U2^ ^ Is once again eliminated by the use of equation (A. 3), and 
the resulting ordinary differential equation for line 2 Is written 
as follows 

( d^u^/dx^ ) + (C^/12 hy ) [16 - 31 U2 + 16 U3 - u^3 

+ (C^/12 h^) t~30u2 + 32ujjy^2"2u2j^^2^ 

• C2 (i + i>2 +(C^/6 hyKdv/dx) | ^ (A. 8) 


Bounding line 3 to NY - 2 do not require the use of any flc- 
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titlous line; however the gynanetry considerations will be used to 
derive the ordinary differential equations. Foi the line NY - 1 
the ordinary differential equation can be written as 


"y > [ - V - 3 * “ “NV - 2 
- 30 j *16 + (Cj/12 ) [ - 30 X 

“nY - 1 “2 NY - 1 “ ^ ^^3 NY - 1^ 


where is the displacement at the fictitious line and 
it is elemlnated by making use of the shear stress boundary 
condition at that face 


yx 


face II 


3v/3x + 3u/3y ■ 0 (A. 10) 


The equation (A. 10) can be expanded as 
'Vf^ ' \y - 1 . dv 

'N ■ ■ KV 

" V*’' ■ “nY - 1 • 2 *’y NY 

f V 

The elimination of u^ from the equation (A. 9) with the 
help of equation (A. 11) results in the following ordinary differ- 
ential equation for the line Uj^^ ^ 

( . i/dx^H(Cj/12 ) [. 3 + « “NY - 2 
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“nY - 1 '*HY ^ ■♦’^Ci/12 h^) [-30 Ujjy . ^ + 

“2 HY - 1 ‘ ^ “3 NY - 1^ " ^2 “ST *^NY - 1 

-(C^/6 hy)(dv/dx)|^ (A.12) 

In the similar manner the ordinary differential equation can 
be written for the line NY, except that now a backward finite dif- 
ference equation has to be tised to write the difference equation 
for the partial derivative with respect to y. For example 

W*u^/dxV(Ci/12 . 3 + 4 j + 

‘ “ot - 1 - “ V + < '^I'll''^ ) [ - 30* 

Uny + 32 U2 jjY - 2 «3 iiY^ - (v *»• w)jjy 

(A.13) 

^)nce again, the displacement u^ on the fictitious line 

f V 

NY is eliminated with the help of equation (A. 11) and it leads 
to following ordinary differential eqttatlon for the line NY, 

(d^UjjY/dx^)+(Cj/12 hy )[ - Uj^ ^ ^ + 4 Uj^ ^ 2 + 

“nY - 1 " ^ “nY^ +(Cj/12h^ )[-30 uj^+32u2j^- 
^ “3 NY^ " ^2 "'•’nY hy)(dv/dx)j„ 


(A. 14) 


1S2 


onoMi: ■ 

or NOR QUMiTY 

Tha ordinary diffarantiaX aquationa for linaa MY + I throu|^ 1, * 

2 MY can be written in the alnilar Bumer. For line 2 MY 4 1 

throng ai^ - 2 MY there is tu> need to take sYvaatry into »»slder'> 

ation to write finite difference equations for tlM partial differ* 

entials with respect to s, since there are points avail^le 

in the grid to write central differ«ice equatiras. 

For the lines - 2 MY 4’ 1 JL^ in addition to 

shearing stress o , the shearing stress a is also incorporated 
xy sy 

into the ordinary differential equations. This is done to eliain- 
ate the displaceaents appearing due to the use of fictitious linM 
while writing difference equatitms for the partial iterivatives 
with respect to z. For exaaple 

^zz ^ jt“ NY + 1 * ® + 3w/9x » 0 (A. 15) 

For the fictitious line NY + 1, the above equation can be 
written as 

fz I 

“ .Z- NY •» 1 ~ 2NY -t- 1 dw 

2 h “ dx NY + 1 

z 


or 


fz 

'* X- NY + 1 


dw. 


“X- 2 NY + 1 ^dx^ 1 - NY + 1 


(A. 16) 


fy 

Sinilarly by making use of the fictitious line u ^_2NY + 1, 

the shearing stress c is reduced to the following equation* 

xy 


m 


S* 255 «-' Mae Bt 
fOtm QUynjfy 


J m ♦ 1 ■ "i- 2 » ♦ J + * S “iH JL- 2MT + 

(4.17) 

Th» ordintry dlff«r«Dtial •^niatloQ for tlio 11m jt- 2 RY <f 1 
is vrittoa boloif, 

< '*"jU 2 R + i/*t*)+<VI2hJ ) [ U , w + 1 - 

“*^-2R.l ♦‘"Jt-2R + 2**“JUJR+J. 
"X- 2 R . 4] +>'i/12 h* )[ - u^_ 4 R ♦ 1 + !•* 

"i- 3 R + 1 * * “X- 2 R + 1 + « »X- R + 1 ■ 

“ R + ®2 IT <’ ♦ "JX- 2 R + 1 <*•“) 

21tata.lt«« of u'5'^ 2 R + 1 1 fw tl» 

•quotlon (4.18) viUi tho tiolp of oquotiono (4.16) and (4.17) 
rotulto In tho followint ordlnory difforontlal oquation, 

< 2 R + i/«i»*)+<Ci/12hJ ) [ -20 «j|. , „ ^ j + 

“ -t- 2 R + 2 + * “X- 2 R + 3 • "X- 2 R ♦ 4> 
(e^/12h* ) t - »X- 4 R + 1 ♦ 2‘ »i- 3 H, + 1 - 

"X- 2 R ♦ 1 “ "i- R - • ®2 4 <» ♦ *)X- 

*<V‘ V*''**X- 2 R + 1 •>,)» 


HY 1 


(dw/dx) jf. iiY + j 


(A. 19) 
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The Mae Mthod is used to cmistruct the ordinary differen- 
tial equations for the reaaining lines. Appropriate difference 
equations are used to express the partial derivati>^. For exsBq>le 
for lines NY 1 through JL* to write the difference equation 
with respect to s, a backward finite difference formula has to be 
t»ed. Since there are not enough gtid points to write ctt&tral 
difference equaticms. 
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A. 2 StrlvfttioQ of Or^Mvy MffovMtiol Iquotlooo for y - Mroo- 
tloBol Xdnoo 

Tho first y-diroctioool lias is fotaod by tiis iatorsoetioa of 
y-s sad y-x coordinsto pisnos (m^ro 7)* Tho riiosr strsssos <m 
thoss plsaos sro utilisod to sliainsts tho fictitious linos. Ibo 
sh«ar stress for tho free surfsM y»s is 

Oyj^ I coordiasts pirns • 0 (A.20) 

or 


9u/3y + du/3x ■ 0 


(A. 21) 


lins 1 


Iteing ths lasginsry displscsasnt 
ths squstion (A. 21) is vrittsn ss 


of ths fictitious 


or 


V -H " ^1 

2 h_ 


fx 


fx 




’i •''w + i** V*''‘**> 


(A.J2) 


Ths plans x-y is ons of ths synstry pisnos mad v - displscs- 
asnts srs syaastric vith rsspset to this plans. Ths fictitious 
displscsasnt is ussd to writs a forward finits diffarsn^ 
squatiem for ths partial dsrivativs with rsspset to x. For ths 
partial dsrivativs with rsspset to s, a Mntral dif formes aqua- 
tion is ussd by using two fictitious displaosasnts md 
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ordinary differential equatimi for line 1 ie erltt«i am 

(d^v^/dy^H(C^/12 ) t 11 Vj** - 2® ^ V + 1 ■*■ 

^^2N2 + l“^3Itt-»- \ ) [- + 16k 

^1^^* - 30 Vj + 16 vj - Vj3 - Cj ^ (u + w)j^ ^^ 23) 


Dm to eywMttry 


fix 

V * V 

1 3 


and 


f2z 

V 18 V 

1 2 


(A.24) 


Substitution of equations (A. 22) and (A. 24) into the «iua- 
tion (A. 23) leads to 

(d^Vj^/dy^-l'fcj^/U h? )[ - 20 + 17 + 1 + 

^2 NZ + 1 ~ ^3 NZ + i^‘*'^Cj_/l2 ) [- 30 Vj^ + 


32 Vj - 2 V, ] - Cj ji 


( du/dy) 


• • 


(u + w)j^ -(11 Cj^/6 hj^)x 


(A. 25) 


In the similar manner by making use of fictitious line and 
the appropriate finite difference equations, the ordinary differ- 
ential equation can be written for lines 2 through NZ - 2. 
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To wrlto tlio ofootiim for Ilao IB - 1, tboro is o of 
flctitlMo liao la tlM I dlroctlon. Lot uo donoto tlw dloploco- 
ant on this tiictttimm lino by • Thio dioplocoMot oill bo 
oliolnotod by arttlnt of OhMr otcooo on tho froo oorfMOt 

^yt I boondliB plooo III * ^ 
loate to 


8»/8o 4 8v/dy ■ 0 


(A. 26) 


oo oxpM^lng tbo oqootion (A. 26) ot lixM HZ, no bovo 


- V 

MZ ^|g 


TT 


- 1 


dw 
dy I bZ 


or 


^NZ 


’•« - 1 - * >■.< 


NZ 


(A. 27) 


In tho sloHor aanI;^r tho dloplocoaont on tho flctitlouo lino 
le - 1** io tenotod by v ^2 _ i • Since 


Bounding plono IV 


loodo to du/2y -f 9v/dx ■ 0 


(A. 28) 


on expanding tho oquotlon (A. 28) ot Use MZ - 1 loodo to 





I 


♦ 


or 
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''nZ - 1 “ ^2 HZ - 1 HZ - 1 (A. 29) 


on Biking use of these two fictitious IIims the ordiniry equation 
for the line HZ • 1 can be written as 

<■*^2 - -»* 
'^NZ-1'*’^'^2NZ-1'‘'^^3NZ-1"^4NZ-i3 ■*“ 

(Cl/12 h^) t - Vjj2 _ 3 + 16 Vjj2 . 2 “ ^ '^NZ - 1 

^NZ “ '^NZ ^ “ ^2 “*■ ”^NZ - 1 (A»30) 

f X fz 

The ficticious displacenencs _ 1 and ellsln- 

ated by substituting equations (A. 28) and (A. 29) into the equation 
(A. 30) and the resulting equation Is 


^NZ - 1^**^ )+(Ci/12h^ ) [ - 20 Vjj2 - i ^2 NZ - 1 
^ NZ - 1 “ "^4 NZ - 1 ^ +(Ci/12h^ ) ^ _ 3 + 


• • 


'^NZ - 2 " "^NZ - 1 ■*■ ^NZ ^ " ^1 dy 


-(11 Ci/6 hjd(du/dy) 


NZ 


^ 1 -(Ci/6 h )0w/dy) 


NZ - 1 


NZ 


(A. 31) 


Using Che method described above, the ordinary differential 


equations for the lines NZ through n - 2 NZ can be written. An 
appropriate choice of fictitious line and the difference equation 
will always lead to a ordinary differential equation. From this 
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•quAtion tb« displAcaant due to f letitioito it oliaiMtod by Mii^ 

tho oliMriiit stroM at tbo boyndary aurf^a. For iiirli, for 
llnaa a - 2 RZ •f 1 throui^ a, aa naad to utlllsa ihaarint atraaa 
<m tbM fraa boiiadi^ plaaa 1. In addltl<m to thla» for IIum 
a - HZ - 1, a - RZ, a - 1 and a aa naad ^dltlonal ^^tiona 
iddch ara provldad by tha ahaarlng atraaa on tha fraa bmmdli^ 
aurfaca on faca III* 

Tha ordinary dlffarantial aquation for tha llna a can ba 
vrlttan aa 

(d*»_/4y*)+(C^/12 hj )t- . aa ♦ * V, . + 

“ V**] *■ - 3 

+ *v^_j + 6v^_l-20»^ + U V,** ]- Cjx 

(i + i>, (A.32) 

Nota tha uaa of backward flnlta dlffaranca aquatlona for 

fx fs 

both tha partial darlvatlvaa. Tha dlaplaeaaanta v. and 

B m 

on tha flctltloua llnaa ara allalnatad froa tha aquatlm (A.32) 
by aaking iwe of tha alMarlng atraaa at tha fraa aurfaca in tha 
following aannar, 

^yx I plana I * ® 


8u/8y + 2v/3x ■ 0 


(A. 33) 
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The equation (A. 33) can be written aa 


or 



- V 


■ - N2 


2 h_ 




SlBllarly 



- HZ hjj)(du/dy 


-(2 h^)(dw/dy) 


(A. 34) 


(A. 35) 


With the help of equations (A. 34) and (A. 35) the ti#o die- 
placeoents which do not belong to grid diaplaceaanta are elinin- 
ated and the new ordinary differential equation for the line m is 


U\/i„h+(C^n2hl >[-»„. ,S2 + ‘ V, . 2SZ + 
” \ M - »> “I 3 ♦ 


(11 Cj^/6 h^Hdu/dy) 


n 


+ (11 hj)(dw/dy) 


n 


(A. 36) 
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A. 3 Dtrlvfttion of Ordinary Mfforontiol Iquotlons for s - Oiroc- 

tionol Linoi 

Am ■hoMn in Fifuro 8, tht llnoo 1 throuih N3CC oro on tho 
crnckod fOM vhilo lima MXC 4> 1 throuilt MX art on tba uneraekad 
faca* By aaklng taia of fletltioua linaa tha ordinary diffarantlal 
oration for llna 1 can ba vrittan at 


(d^w^/da^-KCj^/12 by ) [ 11 + + 

^ *2MX + 1 ** ''3HX + 1 ^ •KCj ^/12 ) Eli - 20 + 

6 W2 + 4 W3 - - Cj (Ii + i)j (A.37) 

tftiara ^ and ara tha diaplacaaanta on tha fletltioua 
llna 1^ and 1^*. Tha ahaar atraaa 


X* 

laadi to 


plana IV 


- 0 


(A. 38) 


8u/9s + dir/dx • 0 


(A. 39) 


Tha ab«‘)va aquation la axpn&dad for tha llna 1 aa followa* 



du 

df 


1 



-f 2 (du/ds) 


1 


(A.40) 


Sladlarly by naklng uaa of tha ahaarlng atraaa on faca v(b) 
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VC can write 

Wi^^ - ^ ^ +2 hy(dv/d*)j j (A.41) 

on cubctltutlng the equetlons (A. 40) and (A.41) Into the eq\ietlon 
(A. 37) Iceda to the folloiflng differential equation for the line 1 


(d^w^/d*^)+(C^/12 hy ) [- 20 + 17 Wjjjj ^ j + 

^’'2NX + l“''3NX+l^ +(Cj/12 h^ ) [- 20 + 

17 W2 4 Wj - w^ ] ■ Cj (u + v)j^ -(11 C^/ 

6 hj^)(du/dz) j -(11 C^/6 h^(dv/dz)L (A.42) 


Siffiilar procedure can be used to write the ordinary differ- 
ential equationa for the linca 1 through NXC. The w diaplacenenta are 
eynsaetric with respect to the plane v(a)« This consideration 

helps in ellnlnating the displacements on fictitious lines. For 

fly 

example by making use of two fictitious lines 
f2v 

^ in y-direction the ordinary differential equation for 

NX — t 

the line NX - 1 can be written as follows 


''nX - 1 ^*** )+(Cj^/12 \ )[- Wjjj^ - 3 ''nX - 2 “ 

^ ''nX - 1 ''nX “ ''nX ^ ‘♦■^Cj/12 hy )[- 

fly f2y 

''nX - 1 ''nX - 1 ■ ^ '^NX - 1 ^®''2NX - 1 


• s 


- w 


3NX 


, 3 - c, (u + v) 
- 1 2 


NX - 1 


(AA3) 
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(A.44) 


and 



" *2 MX - 1 


f X 

Tho displacoMot can hm axpraatad 


'wc ■ "w-1 * ^ K <*’''*»> 


MX 


(A.45) 


With tha halp of tha Miuatlona (A.44) «id (A«45) Uia fictitioua 
fx fly f2v 

diaplacaaanta , ^MX - 1 *'mX - 1 alialnatad fro« 
tha equatimi (A. ) and tha raaulting ordinary diffarantial aqua- 
tion for lina MX - 1 ia 


h* ) [-30 «„ . 1 + 32 vj „ . ^ 

■ * *3 HX - 1 ^ +(C^/12 hjj ) [- »« . 3 + 1* W|OC - J 

■ ''nX - 1 ''mX ^ " ^2 d* ''^NX - 1 “ 

<C^/6 h^) (du/dr) |jj 3 j (A.46) 

SinLlar procadura ia followad to vrita tha ordinary diffar- 
antial aquations for linaa MX 4- 1 throufdi n. 
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iypPEXDIX B 


Application of tbo KocurronM Xolaticm Hithod Bor 

# 

p-Diroctional Linos 


As shown in Figurs 4, tho prssoncs of erock dividos tho 
face V into two facas, naaaly faca V(a) and V(b), raspactivaly. 

Tha faca V(b) is a traction fraa attrfaca. All tha y^liMB starting 
frc» this faca will satisfy tha following boundary condition. 



0 


(I.l) 


On *ha faca V(a) , al? tha prascribad w-displaccsM»ts ara saro. 
HancSp tha y-diractlonal linsi? starting froai this faca tdll 
satisfy. 


V 


V(a) 


- 0 


(B.2) 


Using tha strass - diaplacaaent ralation aquation (B.l) can ba 
raducad to 


( dv/dy) 



I y (du/dai + dw/ds 

(B.3) 


The aouations (B.2) and (B.3) can ba assaablad into a 
vactor fora by following tha saqoonca of ordaring of tha y>dirac* 
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t II 

tiooal lioM. Hm vsetors an V. . ta4 ?, . . Tkm 

IfK 1,E 

bactor V at faca V cma ba vrlttn m foUam , 


l*u 

'i.k 

1 1 


'l.k 


l.u 


<B.4) 


2nl 


♦ II 

wiMra V. and V. raprasants tha udmoim of tha vactor. 

To gat the roulning boundary conditions tha tre.etlon cm^tiMi at 
faca II la used. In the case of tanalla loadit^ we baua an ^ipllad 
o at the faca, aa ahoim In ngura 2. TIaing the atceaa-diaplaca* 
■anta relation, «a cm write 

(4v/dy)| „ 


V 

1 - V 


(du/dx + dw/da) 


II 


(B.5) 


The aquation (B.S) are again aaaablad into a vactor fora 
and are written as 


n,k E 

<«+*> II 


(1 + V) (1 - 2v) 
<1 - V) 


(B.$) 


conplata vector at faca II can now ha written as 


in 


ar * S £; ■ 

fiEfiOQR omyiY 


n»tt 


n,k 


(1.7) 


Vhere, V r«p 

QfU 


its ths pmX of ths VSCtOT. FoUowll^ 

the recurnmcc relstion SMithod, ask eqaetioa slailar to 
tion (3.96) is vrlttvft bslov 




(B.8) 


The evalutstioa of Mtrix and F. is described in Chapter 

n n 

2. For conyenience natrix D and vector F is partitioned «id the 

II II 

equation (B.8) cm be rewrittm as 


n,u 


'n,k 


^11 ®12 ®13 ®14 
^22 ®23 ®24 


V,* 

l,u 

V* 

l.k 

V" 

l.k 

vv 

l,u 


✓ N 


(B.9) 


Where order of various nstrices is 


’n,. 

sod 

Vi 

l,u 

■ *1 ''i.k 

’'.A 

BKl 


{«').! V" 

“n ■ 

t 

^13 

t 

®21 ’ 

1 

*23 

n'. 


BXl 


(M )xl 
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e„ -(«) 

6„ ) 

1 

1 

®22 > 


■' • IBx HXC 


TIm equation (B.9) ia expanded aa. 

t 

t 

’n.k ■ »21 



tt 


®24 ^l»u ■*“ 


The above equetlon could be reerrenged ee 



With the help of (B.IO) all the tnknoims In the vector 
can be detentined. Oncie the couplete vector la known, using 
the recumnee relation aethod dlsplacoaent v and v are obtained 
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•t all th« grl4 points. 

In the abMaco of nonMl strooMo at faea Zt of ^a i^aei* 
aan, tha a<tuatl<m <B.5) nai^ to ba aodifiad by Mtt^i# a a^l 
to zaro. Anothar ^ing whi^ daaarvas attMtiM at this point is 
tha adjustnent of aquati<m (B.9) for tlM situation in idii^ enrvad 
crack fnmt la pmant. In siKh cases tl^ lines are divided in 
two seta. TIm first sat satisfiM the boundary osoditiem givm 
by the aquation (B.l). Total nunbar of th^a IIom is given ^ 
NFREE. Tha second set of liMs satisfy tha equati<» (B.2) at its 
starting boundary ud total msd>ar of sudi linM is given by BFIX. 
Once again tha go^ming equation is slailar to the aquatim (B.9) 
but now the vector can not ba partitioned as ^own in tha 
equation (B.9) . The distrlbutlcni of the known and unknown ale- 
nents of vector is dapendoit tq>on the way it is decided to 
fix particular grid point on face V, idiich in turn depends upon 
the shape of tha curved crack front. To solve Uils new Matrix 
equation, the vector is rearranged so that it could be writ- 
ten as in the aquation (B.9). Now a* is no longer equal to 
NZxNXC, but its new value is equal to NFREE. The diange in the 

arrangenent of V. vector requires a rearrangenent of 0 astrix to 
X n 

maintain the sane original equations. This is achieved by inter- 
changing the coluans in aatrlx. For exaaple, if v^ is inter- 

changed by Vj in vector, than the colusn i in aatrix will 
ba replaced by its J th colusn «dille in the place of coluan J the 
eleoants of column i ^lll be placed. This adjiwtaent gives rise 


to A aodifiod voctor. Tho of this proe^or* 

is that after this rearraogttnt, exactly the 9mm proe«ieini cam 
be followed aa e3q>laliied before for the atrai^t crmt3t. Once the 

f 

■odified vector V im kiK)wn» wltii the help of inforaation for 

each Individual grid point on face v tlw vector la conatructed 

« 

Finally it la uaed to obtain aolutimia for v and v for all tte 
grid points ^ing tiie equaticm (2.91) of the recurr^ce relatl<m 
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APPBMDIX C 

Deriv«tl<m Of lh« Differential Equatlona For 
y>Directional Lime For %eer Loading 

Die y-directi(»al liMt iAmo equations arc directly affected 

by the iriiear loading are 1, 2, .....NZ «id NZ 4 1, NZ 2, ....2 NZ. 

The shear stress at fa^ IV is 
yx 



^yx I IV * ” *r(y) (C»l) 

where t(y) is expressed as 

T(y) ■ o (y/i») (1-y/h) (C.2) 

is the asxinun value of shearing stress at y ■> h/2, and 
h is the seiBi length of the specimen as shown in Figure 1(b). 

Using the stress~displacement relation, the equation (C.l) 
can be written as 

^Ou/3y) 4 Ov/3x)^ ■ - T(y)/C (C.3) 

G is the shear modulus of the material. Using the procedure 
followed to reduce the second Navier - Cauchy equation to an 
ordinary differmtlal equati<m, the modified equation for the 


line 1 is written as 
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+(C^/12 h^ a 11 “ 20 Vj + 6 + 1 " 

^ ^2NZ + 1 " ^3H2 + 1^ -KCj/12 ^*+ 16 

^16v2-V33-C2 + ^ (C.4) 

Dae to eyawtry 



£22 


(C.5) 


The equation (C>3) can be expended for line 1 in the follow~ 
ing aanner. 


or 


fx . 

V 41 
2 h 



T(y) 


* ^NZ + 1 ■'■^2 h^)C(«iu/dy) ^ + (T<y)/G)] (C.6) 

fx fXz f2s 

Elimination of fictitious lines end fr<xn 

equation (C.4) with the help of equations (C.S) end (C.6) leads to 


(d2v^/dy^) (C3/I2 ^ C- 20 + 17 ^ 3 + 

^ ^2NZ + 1 ’ "'m + !]■’■ ^0j_/12 hj ) t - 30 


202 OWGINAL PAGE rS 

JK POOR QUALITY 

32 V2 - 2 V3] - C2 (u + Cr) ^ - ll(Cj/6 )x 

(du/dy) ^ -(11 C^/6 )(T(y)/G) (C.7) 

In a similar manner equations for lines 2 through NZ can be 
modified. A comparison of equation (C.7) with the equation (2.30) 
shows that the two equations are nearly the same except the addi- 
tional term of -(11 Cj^/6 h^)(x(y)/G) which appears due to a non- 
zero prescribed shearing stress a . The modified equations for 

xy 

line 2 thou^ NZ also Include the additional term -(11 C^/6 h^) 
(T(y)/0). Other terms for these equations remain the same as before. 
The equation for line NZ + 1 is written as 

(li ^ j/dy )+(Cj^/12 )[- + j ’'nz + 1 

" ''iJZ + 1 ''nZ + 2 ■ ''nZ + 3 ^ ■'•(Cj/12 )x 
- + 16 Vj - 30 ^ ^ + 16 ^ j - »3^ ^ J 

■ S“df" '“■""^NZ+l • 

Due to symmetry of v displacements 
f2z 

'^NZ + 1 “ ^NZ + 2 

flz 

^NZ + 1 “ ^NZ + 3 (C.9) 


On substituting the equations (C.6) and (C.9) Into the equa- 
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tion (C.8) leads to the following ordinary differential equation 
for the y-directional line NZ -l- 1, 

^ + 1^**^ ^ (^1^^^ h^ ) [ » 30 + 1 '*■ 

^NZ + 2 “ ^ '^NZ + 3^ +(Cj^/12 h^ ) [ 16 Vj - 

^NZ + 1 ‘‘‘ ^2NZ + 1 ~ "^3NZ + 1 ^ 

“ ^2 ^ ''^NZ + 1 '*■ + (Cj^/6h^)(T(y)/G) 

(C.IO) 

Once again the difference between equation. (C.IO) and the 
one used in section 2.2.2 is only the additional term (C^/6 h^)x 
(i(y)/G) , appearing due to the prescription of a non-zero shearing 
stress. The equations for lines NZ + 2 through 2NZ will also have 
this additional term of (Cj^/6 h^) (T(y)/G). 

Introducing the matrix notation, all the ordinary differen- 
tial equations along the y-dlrectlonal lines are expressed in the 
form 

“ir" (’] * Cn) ^ 

mxl mxm mxl mxl mxl 


(C.ll) 
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APPENDIX D 


Three Dimensional J-Integral 


The definition of J-lntegral presented in this appendix is 
taken from reference (37). Consider the section of three dimen- 
sional crack frcmt shown in Figure 47. The surfaces, S^, $2 and 
surrounds the crack front. The component of the general 


"r 

J-integral vector J is given as 
J “Jr 

X ''1 


(Mdy - I.J Uj ,, d.j) 




(o. U )/ dS 
Iz i,x 


(D.l) 

A,* ,X 

T2 + S3 

where i « 1, 2, 3 and P is a circular path as shown in Figure 


47. W is strain energy density, is stress tensor and u^ is 

referred to as displacement vector. 

In two dimensional plane stress problems becomes zero 

and the equation for becomes the same as given by Rice (27). In 

the case of plane strain problem o is not zero but its variation 

zz 

through the thickness is zero. Therefore, the surface integral 
term drops from the equation (D.l) and once again the equation for 
becomes the same as given by Rice (27). 
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APPENDIX E 

Description And Listing Of The Coq>uter Programs 

In this appendix a description of the computer programs used 
in the present study is given. A complete listing of all these 
programs including subprograam is attached at the back of the 
thesis. In writing the various programs strict modxilar approach 
was followed. This has helped in writing many subprograms which 
are commonly shared by several main programs. 

The program was divided into five main progranm. This was 
done to divide one complete run into several small runs. This was 
done to safeguard against the failure of coii 4 >uter aystem « The 
schematic representation of the computer program MAINl is shown 
in Figure 48. This program is enq>loyed to evaluate the various 
matrices needed for the X'dlrectlonal calculations. At the end of 
the computer program run, all these matrices are stored on the 
virtual disc storage. Computer programs KAIN2 and MAIN3 %diose 
schematic representation are given in Figures 49>*S0, are also used 
to calculate different matrices used for y and s-directlonal cal- 
culations. The outputs from these programs are also saved on the 
virtual disc storage. 

A dynamic storage allocation scheme was used for all the 
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programs. Consequently whenever there was a change in the nun^er 
of grid lines, only the DIMENSION statements In the main programs, 
need to be modified. Due to the serious round cff error problem, 
all the calculations for the x-dlrection were performed in 
double-double precision (128 bit word size). For the other two 
directions the double precision was sufficient. 

The computer program MAINA Is used to obtain the elastic 
solution for different cases. The subroutines UDISP, CCV, and 
WDISP perform the calculation for x, y and z-dlrectional lines, 
respectively. In the case of tensile loading subroutine CCV was 
replaced by subroutine CVDISP. The subroutine OUTPUT is used to 
print the results for the displacements and their derivatives. 

The normal and shear stresses are calculated with the help of sub- 
routines STRESS, SHEAXY, SHEAYZ. The effective stress at each 
grid point can be computed by using subroutine EFTRES. The output 
from the program MAIN4 was also stored on the virtual disc. 

Four virtual tapes are needed to run the program MAIN4. The 
stored matrices which were originally calculated by programs MAINl, 
MAIN2 and MAIN3 are read into three sequential virtual tapes. The 
fourth virtual tape is needed to write the output of program 
MAIN4 on the virtual disc. 

Program MAINS is used for the elastic-plastic stress analy- 
sis. Its schematic representation is shown in Figure 52. There 
are many subroutines which are shared by this program with the 
program MAIN4. Due to the addition of plastic strain terms many 
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old subroutines were modified. Note the presence of new subpro- 
grams PUDISP, PCCV and PWDISP, used to perform calculations for 
the X, y, and z-dlrectlonal lines, respectively. This program uses 
five virtual tapes. Three tapes are utilized to read the different 
matr.ces generated by the programs MAINl, MA1N2, and MAIN3, The 
remaining other two virtual tapes are used to read and write the 
solution, at the beginning and at the end of the increment. 

This format of programming provided considerable flexibility 
and a very good protection iu case of computer systems failure 
during the execution of the program. In case of a system failure, 
only the solution for the current increment is lost, which can be 
recomputed by using the stored solution for the previous increment. 


fAGg IS 
QUALrrv 


PROGRAH HAINI 
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